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1. Introduction. The theory developed in the present paper applies to any elastic 
body (in general anisotropic) possessing a positive-definite strain-energy function, 
quadratic in the components of stress. This function provides us with a metric in a 
function space in which the point or vector represents a state of stress. The geometry 
of the function space follows Euclidean analogies closely, and is powerful in suggesting 
methods of approximation. The aim is to obtain approximate solutions of elastic 
boundary value problems with errors which are calculable, the error being measured 
in terms of distance in function space, or, equivalently, in terms of strain-energy. 

After dealing with notation and basic concepts in Sec. 2, we discuss vectors in 
function space in Sec. 3. Sections 4 and 5 are intended to introduce the reader to the 
ideas which'lie behind the general method; only a first approximation is discussed, 
and only the simplest types of boundary conditions. 

In Sec. 6 more general types of boundary conditions are introduced, and higher ap- 
proximations under these boundary conditions are treated in Secs. 7-11. Throughout, 
the basic plan is to locate the solution of the elastic boundary-value problem (con- 
sidered as a point in function space) on a hypercircle of determinable center and 
radius. As a practical test of the method, it is used in Sec. 12 to obtain approximate 
solutions for the torsion of a prism of square cross section. In Sec. 13 reference is made 
to other work, and some known results are strengthened by use of the present method. 

2. Notation and basic concepts. Latin suffixes take the range of values 1, 2, 3 and 
the summation convention operates on repeated suffixes. The coordinates x; are rec- 
tangular cartesians, and differentiation with respect to a coordinate is indicated by a 
comma (F,;=0F/0x;). 

By a state of an elastic body we understand a set of six stress components Ej; 
(£;;=E,:), given as functions of the coordinates throughout the body. For simplicity, 
we shall consider only functions which are continuous in the body and on its surface, 
and possess continuous first and second order derivatives in the body. It will be obvi- 
ous that these conditions may be weakened, in the sense that these conditions hold 
throughout each of a finite number of parts into which the body is divided, with 
suitable continuity conditions across the surfaces which separate the parts. 

The body is assumed to possess a strain-energy function 


W = 865 jmalijEmn; 
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where the constants Cijmn satisfy 
Cijmn = Cjimn = Cijnm = Cmnijs (2.2) 


and are such that the form W is positive-definite (i.e. it is positive unless all the stress 
components vanish). 

We now introduce the strain components e;; in an unconventional way. Instead of 
defining strain in terms of displacement, we define it by the generalized Hooke’s law 


Qj = Citmakbun (2.3) 


Thus (2.1) may be written 
W = 46; ;E;;. (2.4) 


It follows from (2.2) that for any two states, E;; and Ej;, we have the reciprocity 
relation 


/ In A c 
€:;Ei; = €;;E3;. (2.5) 

The usual equations of compatibility read 
Cij,mn + €mn,ij = Cim,in + Cin,im- (2.6) 


By (2.3), these can be translated into conditions on the stress components. In general, 
they will mot be satisfied by an arbitrary state. But if they are satisfied, then the 
partial differential equations 


Ui; + Uj,i = 2e:;, (2.7) 


have solutions “;, unique to within an infinitesimal rigid body displacement,’ and u; 
is the displacement for the state F;;. We are to remember that an arbitrary state has 
in general no corresponding displacement. 

Should we wish to construct a state satisfying the equations of compatibility, 
the simplest plan is to choose a set of displacements, use (2.7) to obtain e;;, and then 
solve (2.3) to obtain E;;. The compatibility equations (2.6) will be automatically satis- 
fied on account of (2.7). 

Throughout this paper, we shall suppose that body forces are absent.* Thus the 
equations of equilibrium read 


Ei;,; = 0. 2.8) 


In general these equations will not be satisfied by an arbitrary state. 
We shall denote by n; the unit vector normal to the surface of the body, pointing 
outward. Then the stress across the surface is 


T; = E; jj. (2.9) 


9 


The following notation for “inner products” will be found convenient. Since it will 
be obvious whether an integration throughout the body or an integration over its sur- 
face is implied, one type of notation serves for both types of inner product. We shall 
write 





* The extension to the case where body forces are present will be made in a later paper. 
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(e . E’) = f e;;E;,dv, 
(2.10) 


(u- 7”) = f urias, 


the former integral being taken throughout the body, and the latter over its surface. 
The interpretation of similar expressions will be obvious to the reader. 

In attempting to solve a problem in elasticity, we seek a state which satisfies 

(i) the equations of compatibility (2.6), 
(it) the equations of equilibrium (2.8), 
(iii) the assigned boundary conditions, which will generally be T; assigned, or u; 
assigned, or 7; assigned in part and u; assigned in part. 

We approach the solution by considering states in which one or more of these condi- 
tions are relaxed. 

3. Geometry in the function spaces of states. Since 
the power of the present method lies largely in the stimu- B 
lation of geometrical intuition, it is well to pictorialize 
from the start. We shall develop some simple properties 
of the function space of states. S"=$+s! 

The unstressed state, E;; = 0, is represented by the 
origin O (Fig. 1). Any other state is represented by a point 
such as A. Just as we describe the position of a point in 
ordinary space by a position-vector drawn from the 
origin, so we describe the point or state A by the vector o 
OA, or more compactly by a single letter (S) in heavy Fic. 1. 
type. Thus the symbol S represents a state of stress (six 
functions of position throughout the body). 

We add vectors in an obvious way. If we have two states, 






S with stress components £;;, and 
. +f 
S’ with stress components £;;, 


we define the sum S’’=S+S’ to be the state with stress components 
” , 
E;; = Ei; + E;;. 
The definition of S—S’ is obvious. 
To multiply a vector by a scalar constant k, we write S’=kS, and define S’ by 
the equations E};=kE,;. Note that E;; are functions of position in the body, but is 
a constant. The distributive laws are satisfied, so that 


(k+kR)S=kS+HS, KS+S) =kS+4 KS’ 


So far we have mentioned only vectors drawn out from the origin O, i.e. bound 
vectors. But the idea of a free vector, familiar in ordinary space, can be used with ad- 
vantage. Apart from the geometrical representation, there is no real distinction be- 
tween a free vector and a vector bound to O: each corresponds to a given set of stress 
components. But in the geometrical representation, a free vector S starts from any 
initial state A and proceeds to that state B for which the stress components exceed 
those of A by the stress components corresponding to S. In fact, a free vector may be 
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regarded as a transition from one state to another, and we make it a bound vector 
if we specify the initial state. 

The algebraic procedures apply to free vectors. The addition of vectors follows 
the familiar parallelogram law (Fig. 1). 

It should be amply clear from the context or from the appropriate diagram 
whether a free or bound vector is meant. In any case of doubt, a vector should be 
interpreted as a vector drawn from the origin O. 

Fig. 1 will also draw the reader’s attention to the geometrical relation between the 
vectors S and RS. If k were negative, there would be a reversal of sense. 

We define tke length or magnitude S of a vector S by 


S? = f 2Wdv = f é:;E;;dv = (e-E), S20, (3.1) 


the integrals being taken throughout the body. Here e;; is the strain corresponding to 
the stress E;;, by (2.3). Since the strain-energy W is positive-definite, the length of 
any vector is real, and is zero only for the unstressed state, E;;=0. Sometimes it will 
be convenient to denote the length of a vector S by |S|. A vector S is a umit vector 


if S=1. 
The distance between the extremities of two vectors, S and S’, drawn from O, 


is defined to be the length of the vector S—S’. Thus this distance is given by 
|\S—s’|?= f (ei; — e:;)(Ei; — Ej;)de. (3.2) 


It vanishes if, and only if, S=S’, i.e. if the two states are identical. 
The scalar product plays an important role. We define the scalar product of two 
vectors as the inner product of the corresponding states: 


S-S’ =f eiskistn = (c-E’). 3-3) 


It is obvious, from the reciprocity relation (2.5), that the scalar product has the com- 
mutative property 
S:S’ = S’-S. (3.4) 
Two vectors are said to be perpendicular or orthogonal if their scalar product vanishes. 
Let us explore the physical meaning of the scalar product. Let us suppose that S 
satisfies the equations of compatibility, so that it arises from a displacement, and that 
S’ satisfies the equations of equilibrium. Then 


S-S’ = (e-E’) =f eisBisd = fm. sBhido 
=f uBinds — fF uBisds, (3.5) 


by Green’s theorem. The last integral vanishes, since S’ satisfies the equations of 


equilibrium (2.8). Hence 


S:S’ = f u;sT {dS = (u-T’); (3.6) 
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in words, the scalar product S-S’ equals the work done in the displacement of S by the 
surface stress of S’. If the vectors are orthogonal, this work is zero. 
As an example, consider a prismatic bar. Let the states S and S’ correspond to 
simple tension and flexure, respectively. Then S is orthogonal to S’. 
It should be noted that the scalar product can be interpreted in terms of work 
only if one of the states satisfies the equations of equilibrium and the other the equa- 


tions of compatibility. 
We come now to a fundamental inequality. From the positive-definite character 
of strain-energy, it follows that if S and S’ are any two arbitrary states, and k an 


arbitrary real number, then 
fic = ke; ;)(Ei; = kE;;)dv Pa 0, (3.7) 


the integral being taken throughout the body. This reduces to 
S’? — 2kS-S’ + k*S? 2 0, (3.8) 
and since this holds for all real values of k, we deduce the Schwarzian inequality 
|S-S’| < SS’. (3.9) 
If the sign of equality holds in (3.9), then it is easy to see that the states S and S’ 
must be connected by a relation S’‘=KS, where K is some real constant; in other 


words, the vectors S and S’ have the same direction, or opposite directions. 
We can now define the angle 0 between two vectors S and S’ by 


S-S’ 


P. 6S. 3.10 
= a (3.10) 


0 


IIA 





cos § = 


By virtue of (3.9), the angle so defined is always real, and is of course equal to 37 
when the vectors are orthogonal. For the angle between two unit vectors, I and I’, 
we have 


cos@ =I-I’. (3.11) 


The whole of Euclidean geometry holds in any linear subspace of our function 
space based on a finite number of vectors. In particular, we have the theorem of 
Pythagoras, the theorem that any side of a triangle is less than the sum of the other 
two sides, and the theorem that the greater side of a triangle is opposite the greater 
angle. All this is well known, but the abstract and general form usually given to the 
theory of function space may easily obscure the intuitive simplicity of the approach. 

Scalar multiplication is distributive, as is easily seen, so that 


S:-(S2 + S3) = §,-S. + S;-S3. (3.12) 


As an exercise, let us prove the theorem of Pythagoras. We have, for any two vec- 
tors, S and S’, 
(S — S’)-(S — S’) 
= S? + §”? — 28-S’ 
= §? + S’2 — 2SS’ cos 6, (3.13) 


1s — 8’? 
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where @ is the angle between S and S’. Putting 6 = 37, we get the theorem of Pythag- 
oras. 

4. First approximation with surface stress given. Let us suppose that the surface 
stress T; is assigned, satisfying the conditions of statical equilibrium, so that a solution 
S exists. We may refer to the solution S as the natural state, to distinguish it from the 
artificial states which we shall introduce. 

Let S* be a state which satisfies the equations of equilibrium and the boundary 
conditions, but not the equations of compatibility. (If the equations of compatibility 
were satisfied, we would have S*=S, the natural state.) Let us evaluate the scalar 
product S-S*. We find 


S-S* = (e-E*) = f eisBise = f wsBido 
= f “;T;dS = f u;T dS = (e-E) = S?; (4.1) 


in carrying this out, we have used the fact that S* satisfies the equations of equilib- 
rium, and also the fact that it satisfies the boundary conditions, so that 7* = 7;. Since 
S?=§S-S, we may write (4.1) in the form 

S-(S — S*) = 0. (4.2) 
This tells us that the vectors S and S —S* are orthogonal; this means that the extremity 
of S is located on a hypersphere H having the vector S* for diameter; in fact, (4.2) is the 
equation of this hypersphere, S being regarded as a current vector. The center of the 
hypersphere is at 3S*, and its radius is 3.S* (Fig. 2). The equation (4.2) may also be 


written in the equivalent form 


|S — 38*| = 35%, (4.3) 
s" 
ice 
° 
Fic. 3. 





which shows up the center and the radius. 

We now take another state S” which satisfies the equations of compatibility, but 
not in general either the equations of equilibrium or the boundary conditions. In fact, 
we can build such a vector S” by choosing an arbitrary set of displacements through- 
out the body. Let us investigate the scalar product S-S”, S being the natural state as 
before. It is easy to see that 


S-S" = (c"-E) = (u"-T). (4.4) 
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Although S is unknown, we do know 7; from the assigned boundary conditions, and 
so the expression (u”- 7) is calculable. If @ is the angle between S and S”, we have 


S cos @ = §-S"/S” = (u”-T)/S’. (4.5) 


The right hand side is a definite calculable number, and S cos @ is the orthogonal 
projection of S on §”. It follows that the extremity of the natural vector S lies on a 
definite hyperplane P which is orthogonal to S” (Fig. 3). 

Now 


S*.S" = (e"-E*) = (u"-T*) = (u"-T7). (4.6) 


Comparison with (4.4) reveals a rather remarkable fact: the extremity of the vector S* 
lies in the hyperplane P. 
On account of (4.6), we can write (4.4) in the form 


S:S’ = S*-S’; (4.7) 


thinking of S as a current vector, we may regard 
this as the equation of the hyperplane P. 

We have now located the extremity of the 
natural vector S on the hypersphere (4.2) and the 
hyperplane (4.7); we have therefore located it on 
the hypercircle T which is the intersection of the 
hypersphere and the hyperplane. The situation is 
therefore as shown schematically in Fig. 4. The 
hypercircle IT passes through the extremity of S* 
and also through 2”, the foot of the perpendicular 
dropped from O on P. We have 





s*.” 
sm 





>" _ Ss’ (4.8) 


If I” is a unit vector codirectional with S” (and therefore also satisfies the equations 


of compatibility), we have 

I’ = S’/S’, 
and (4.8) can be written in the form 

>” = I’(S*-1’). (4.9) 
The points 2” and S* are at the ends of a diameter of the hypercircle I, because 

(S — 2’)-(S — S*) = 0 
as follows from (4.2) and (4.7). The center C of T is therefore at 
C = 3(S* + 2”) = 3[S* + 1(S*-1’)], (4.10) 

and the radius R of T is given by 


R = 4(S*? — BYU? = B[s? — (S*-1")2}2 (4.11) 
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Of all points on T’, the endpoint of 2” is nearest to O and the endpoint of S* most 
distant from O. Thus, 2”? $.S? < S*? or 
(3° -F)* s S* = S™. (4.12) 


This relation gives upper and lower bounds for the strain energy associated with the 
state S. 
We may also write (4.11) in the form 


R = 35S* sin 0, (4.13) 


where @ is the angle between S* and S”. 

Let us now consider the question of getting an approximation to the solution S 
of the given boundary value problem. Various definitions of a “good” approximation 
might be given. We might base a definition on the maximum deviation of displace- 
ment or stress from the correct value. We shall, however, use our geometrical picture 
as the basis of a test for goodness, and say that an approximate solution § is good if the 
distance (in function space) between S and the true or natural solution S is small. We 
shall define the error e of any approximate solution by 


«e=|S-S|. 


Returning to the discussion of the hypercircle , we ask whether our knowledge 
based on S* and S” enables us to give an approximate solution. Suppose we accepted 
S* as an approximate solution, how great an error would we commit? It is easy to 
set an upper bound to this error, because S lies on I’, and the point on I’ most distant 
from S* is the point 2”, which is diametrically opposed and at a distance 2R. Thus we 
may state that the error of S* does not exceed 2R, where R is given by (4.11). 

Similarly, if 2”, as in (4.8) or (4.9), is taken as an approximate solution, the error 
does not exceed 2R. 

We can however do better. Jf we take C, the center of the hypercircle T, as approxi- 
mate solution, the error is precisely R. Thus, if we change from S* to C, we exchange an 
unknown error which may be as great as 2R for a certain error which is only R. 

It should be pointed out that the approximate solution C will, in general, not sat- 
isfy the equations of compatibility, nor the equations of equilibrium, nor the boundary 
conditions. 

The hypercircle I’, formed for states S* and S”, gives us a yardstick by which 
we can assess any proposed approximate solution §. If D, is the least distance of § 
from I and D, the greatest distance, then the error e of S satisfies 


Di Se = Dz. (4.14) 
It is easy to prove by elementary trigonometry that Dj and D3 are given by 
(S — C)?+ R? F 2R [(S — C)? — (1"-(S — C))?]#/2. (4.15) 


5. First approximation with surface displacement given. In the preceding section 
we have considered the case where the surface stress 7; is given. We shall now take up 
the other type of problem in which the displacement u, is given on the bounding sur- 
face of the body. The argument runs closely parallel to that of Sec. 4. 

First we choose a state S* which satisfies the equations of compatibility and the 
boundary conditions. This is a much easier task than we had in Sec. 4, because 
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all we have to do is to take a system of displacements which assume assigned values 
on the bounding surface, and satisfy the general conditions regarding continuity of 
derivatives. 

Next we choose a state S” which satisfies the equations of equilibrium, but not in 
general the equations of compatibility, so that satisfaction of the boundary conditions 
is meaningless, since displacement does not exist. Let I” be the corresponding unit 


vector. 
We have 
S-S* = (e*-E) = (u*-T) = (u-T) = S?. (5.1) 
Accordingly, the endpoint of S lies on the hypersphere 
S-(S — S*) = 0. (5.2) 
Also, 
S:S’ = (e-E”) = (u-T”) (5.3) 
which is calculable, and 
S*-.S” = (c*-E”) = (u*-T’) = (u-T”). (5.4) 
Thus, 
S-S’ = S*-S’. (5.5) 


Let us compare the equations (5.2) and (5.5) with (4.2) and (4.7), respectively. 
They are formally the same equations, although S* and S” now have new meanings. 
Hence, just as we obtained (4.10), (4.11) and (4.12), we are now led to the result: 
when the surface displacement is given, the natural vector S has its endpoint on a hyper- 
circle with the center at 


C = }[S* + 1'(S*-1’)], (5.6) 
and the radius R given by 
R = 3[S*? — (S*-1")?]!/2; (5.7) 
the strain energy associated with the state S is bounded in accordance with 
@-P)* 4S = S™. (5.8) 


6. Boundary conditions. Two types of boundary conditions may be regarded as 
fundamental: (i) stress assigned, (ii) displacement assigned. As we have seen in Secs. 
4 and 5, there is remarkable duality between these two types of boundary condition. 
This duality extends to more general types of boundary condition, which will now be 
described before we proceed to generalize the ideas developed in the preceding sec- 
tions. 

We consider boundary conditions which involve the following elements at points 
of the bounding surface: 

normal stress T (n); tangential stress T (4); 

normal displacement un); tangential displacement 4 (1). 

The normal elements are of course scalars and the tangential elements vectors in the 


tangent plane. 








W. PRAGER AND J. L. SYNGE [Vol. V, No. 3 


to 
wn 
—) 


We shall classify points on the bounding surface according to the conditions as- 
signed at them. Thus, a point will be said to be of the class [un), T (| if the normal 
displacement and the tangential stress are assigned at the point. Cases of vanishing 
elements are important, and for them we shall use the following type of notation: 
[4cm) =0, T)] means that the normal displacement and the tangential stress are as- 
signed, and that the normal displacement is zero. 

Two types of boundary condition will now be defined. These may be regarded as 
generalizations of the simple types of boundary conditions in which stress is assigned 
all over the surface, or displacement is assigned all over the surface. 

STRESS BOUNDARY CONDITIONS (or briefly SBC): The bounding surface may be 
divided into regions such that all the points in any region belong to one of the follow- 
ing classes: 


[T n) Tw], [T in), u(t) = 0], [Tw), Un) = 0], [1¢.(n) = 0, ut) = QO}. (6.1) 


Obviously, if the complete surface stress T; is assigned all over the surface, all 
points belong to the first class, and so this is a particular case of SBC. Note that in 
SBC any displacement elements which may be assigned are zero. 

DISPLACEMENT BoUNDARY CONDITIONS (or briefly DBC): The bounding surface 
may be divided into regions such that all the points in any region belong to one of the 
following classes: 

[u n)y Ub »], [ tn); Ta) = 0], [¢cy, Ton) = 0], [T ¢n) = 0, TQ) = 0}. (6.2) 


Note that the case where displacement is assigned all over the surface is a particu- 
lar case of DBC, and that in DBC any stress elements which may be assigned are zero. 

As an illustration, consider a horizontal beam carrying a load on its smooth upper 
surface and supported on smooth unyielding supports. The points of the surface may 
be classified as follows: 

upper surface under load [7T(n), Ti) =0] 

surface in contact with supports [T(.=0, un) =0] 

rest of surface [Tn =0, Ty =0] 

The boundary conditions are evidently SBC, only the first and third classes of (6.1) 
occurring. 

As a second illustration, consider the classical Saint Venant torsion problem for a 
beam of arbitrary section. The classification is as follows: 

ends of beam [uc1), Tin) =0] 

sides of beam [T cn) =0, Ten =0]| 
Here we have DBC, only the third and fourth classes of (6.2) occurring. 

7. Associated and complementary differential equations and states. When we 
have to deal with SBC, we have reason to regard the equations of equilibrium, 
rather than the equations of compatibility,,as particularly associated with the 
boundary conditions. On the other hand, when we have to deal with DBC, we regard 
the equations of compatibility as particularly associated with the boundary conditions. 
The real reason for these associations lies in the theory which follows, but we may 
note the following facts which will help in remembering which differential equations 
are associated with which boundary conditions. 

In the simplest type of SBC (E;jn; assigned) and in the equations of equilibrium 
(E;;,;=0), the stress components are involved in a very simple way. On the other 
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hand, the assignment of DBC would be meaningless unless the equations of compati- 
bility were satisfied, because displacement would not exist. 

We shall therefore speak of associated and complementary differential equations in 
accordance with the following table: 














TABLE I. 
Type of boundary Associated differ- | Complementary 
condition | ential equations | differential equations 
| a a ee ee 
— Equilibrium | Compatibility 
| Stress (SBC) | (2.8) (2.6) 
Displacement (DBC) “aa ~~ 








This terminology will permit us to discuss both types of boundary condition with a 
single argument, which bifurcates only for details. 

We note that the solution of an elastic problem involves the satisfaction of 

(1) the associated differential equations, 

(2) the complementary differential equations, 

(3) the boundary conditions. 

The following table sets forth schematically notation and terms to be used: 

TABLE II. 








= | | 














| Stress Boundary Conditions | Displacement Boundary Conditions 
citi Differential | (SBC) | (DBC) 
ine uene Name Equations |— so [ 
eT Satisfied Diff. Equations | Boundary Condns. | Diff. Equations) Boundary Condns, 
Satisfied | Satisfied Satisfied Satisfied 
Associated Equilibrium Hee rol | Compatibility | joes wal 
S Natural state and and | BA Cm> H(e) ™ and | Mn). 2) = 


cabins TQ), =0 | ee » Tim) =0 
Complementary; Compatibility fhe ge yh : ] 0} | Equilibrium | iT o T% 1 0] 
| | (n) =Y, Ut) = | (n) =Y, a= 








| Ti’) =Tm) where 
T (n) assigned for S; 


u(%) =U(n) where 


Completely . . 
ompletely u(n) assigned for S; 























Ti) =0 where 


(¢g=1,°°* ,%) states | uc,) =O where | 
T(t) =0 for S. 


uct) =0 for S. 


S* associated Associate Equilibrium | Compatibility 
aor ssociated | qu T() =T«) where pi ty uc?) =u) where 
| T(t) assigned for S. u(t) assigned for S. 
‘ree << ~~ ae ae re ’ 
(n) =O whe | =0 
, Homogeneous | | | 7 =) wie = Tn) | ¥O) where sgn) 
V8 ‘ . | oe | assigned for S; - — assigned for S; 
, associated Associated | Equilibrium | 7 - | Compatibility ’ 
(p=1,°°**,m T «@) =0 where T(¢) ui) =O where u (4) 
states | | 4 4 
| assigned for S. assigned for S. 
| u(,,) =O where | T(,.) =0 where 
i Complementary}, " ce u(n) =O for S; ae Tn) =0 for S; 
¢ = ¥/Complementary| Compatibility i : Equilibrium | () 





Orthonormal ————_- 














i, homogeneous -| As for r aS area —? 
(p=1,°**,m)| associated — ' 
states 
7” Orthonormal | — 7 
, . ) complementary| *-| As for S, a - : anemone memenenteied 
7 coo lg imal? 


states 
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We note that §, the natural state, is the solution of the elastic problem with which 
we are concerned. It is to be regarded as unknown. All the other states are artificial 
states which we build up with a view to getting an approximation to §; they are 
to be regarded as known. 

In any particular problem, the boundary conditions will be either SBC or DBC. 
Hence, in any particular problem, we shall have occasion to consult the two columns 
under SBC, or the two columns under DBC, according to the type of boundary con- 
dition. 

The table is self-explanatory, but it may assist the reader if we interpret it in the 
case where the stress 7; is given all over the surface. Here we have SBC, and the 
following interpretation: 

S* is any state which satisfies the equations of equilibrium (2.8), and has at 

each point of the surface the assigned value of T;. 

» form a set of states, each satisfying the equations of equilibrium and homo- 
geneous boundary conditions T;=0. 
S/’ form a set of states, each satisfying the equations of compatibility, but not 
subjected to any boundary conditions. 
» form an orthonormal set of states, each being a linear combination of the 
states S, ; each therefore satisfies the equations of equilibrium and makes 
T;=0 on the surface. 

I/’ form an orthonormal set of states, each being a linear combination of the 

states S/’ ; each therefore satisfies the equations of compatibility. 

The procedure for forming an orthonormal set by linear combination of linearly 
independent states is a well known procedure. The general conditions for orthonor- 


mality read 
In-l=5p (p,r =1,---,m) 


47 Ad ‘Oe 
I, I, = 5¢s (q,s =1,---,m). em 


These conditions imply that each vector is of unit magnitude, and that all vectors 
in each set are mutually orthogonal. 
We shall always assume that the states S/ are chosen linearly independent; other- 
wise the number of orthonormal vectors I, would be less than the number of vectors 
,. The same remarks apply to the states S/’. 
The following statements are easily verified: 
I. Any linear combination of homogeneous associated states is itself a homogeneous 


associated state. 
I]. If any linear combination of homogeneous associated states is added to a completely 


associated state, the resulting state 
S* + >0 4,8, (7.2) 
p=1 


is a completely associated state. 
III. Any linear combination of complementary states 
de a, S, (7.3) 
q=1 


is itself a complementary state. 
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8. The associated and complementary hyperplanes. If, in (7.2), we assign arbi- 
trary values to the constants a/ , we get a linear subspace or hyperplane of m dimen- 
sions. We shall call this the associated hyperplane, and denote it by L},. Similarly, 
we define the complementary hyperplane by (7.3), and denote it by L,’. Fig. 5 shows 
these hyperplanes for the case m =1, n=1, when, of course, they are straight lines in 


function space. 








Fic. 5. Fic. 6. 


The associated hyperplane and the complementary hyperplane play important 
roles. If they possess a point of intersection, then the state corresponding to that 
intersection is the natural state S; for this state satisfies the associated differential 
equations, the complementary differential equations, and all the boundary condi- 
tions. The object of our approximations is to draw these hyperplanes together, and 
to locate their closest points, so that we may get a good approximation to S. 

We shall now prove the following: 

IV. The associated hyperplane LX, and the complementary hyperplane L,' are mu- 
tually orthogonal. 

The hyperplanes are orthogonal if every vector lying wholly in one hyperplane is 
orthogonal to every vector lying wholly in the other. The general vector lying wholly 
in L,,* is 


m 


ty! 
S' = > a,l), 
p=1 
v? 
and the general vector lying wholly in L, is 
n 
cas Whett 
S’ = ) a l,. 
q=1 
The necessary and sufficient condition for the orthogonality of the hyperplanes is 


s’.S’ = 0. (8.1) 


To prove that this relation holds, we must consider SBC and DBC separately. For 
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SBC we have, since S’ satisfies the equations of equilibrium and S” the equations of 
aon oe 0 os 
S’.S’ = (e’.E’) = ff eiziae = fo, Bisdo 


| a Phent 
= f u; T jd - fu E;;,;dv 

eT sal . ae - : 
= f WisTinas +f ur) 1 (dS, (8.2) 


the last integrand being the scalar product of vectors of displacement and stress lying 
in the tangent plane of the bounding surface. Let us now use the columns under 
SBC in Table II, and consider the integrals in the last line of (8.2) taken over the 
several regions into which the surface is divided. We see that 


compatibility, 


on [T <n)» ig | we have Ths = 0, . = 0, 
on [T en), ue) = 0] we have Too = 0, uu) = 0, 
on [Tc hia) = 0] we have Tis = (), ad = (), 
on [1¢(n) = 0, #0) = 0] we have ‘ = (0, u 1) = 0. 


Hence the integrals in the last line of (8.2) vanish, and (8.1) is proved; this estab- 
lishes the orthogonality of the two hyperplanes. The proof for DBC follows the same 
lines. 

From (8.2) we have 


>. se 


I,-I, = 0 (p =1,2,---,m;q=1,2,---,m). (8.3) 
9. Some basic relations. We shall now prove a number of results. 
\V. For the natural state S and any completely associated state S*, the following rela- 
tion holds: 
s-(S — S*) = 0. (9.1) 
Hence the extremity of S lies on a hypersphere having the vector S* for diameter; the center 


is at 48* and the radius is 4S*. 


To prove this, we separate the cases SBC and DBC. For SBC we have 


S:S = (4.T) = f HT ord + f mT dS, 


S-S* = (u.T*) = f wmTdas + f moTdas. 


By Table II, 


sa = an a = 

on [T «ns T «| we have Tm) = Ta), Tay = Tay 
on Uk 

on [Tin), uy = 0] we have Ta) = Ta), Ua) = 9, 
ms ait , 

on [Te), Cin) = 0] we have Ta) = Ta, 4) = 0, 


on [1¢¢n) = 0, 4a) = 0] we have um) = 0, ua) = 0. 
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On using these values in (9.2), we see that S-S =§S-S*, and so (9.1) is proved for SBC. 
The proof for DBC is similar. Since (9.1) may be written 


(S — 3S*)? = 15*2, (9.3) 


the truth of the statement about the hypersphere is obvious. 
VI. The natural state S is orthogonal to any homogeneous associated state S’: 


S-S’ = 0. (9.4) 
Hence S is orthogonal to the associated hyperplane L%, and 
SI,=0 (p=1,2,--+,m). (9.5) 


Since any homogeneous associated state may be expressed as the difference be- 
tween two completely associated states (cf. Table II), (9.4) follows if we write down 
relations of the form (9.1) for the two completely associated states and subtract one 


from the other. 
VII. The difference between the natural state S and any completely associated state 
S* is orthogonal to any complementary state S": 


(S — S*)-S’ = 0. (9.6) 
Hence this difference is orthogonal to the complementary hyperplane L,!’ and 
ois. S*-I, (q = ee ee) (9.7) 


For SBC we have 


S-S” = (e’-E) = (u’-T) = f tiny T (nydS + f wisTwas, 
(9.8) 


S*.S” = (e”-E*) = (u"-T*) = f ol ott + f miTYoas. 


On reference to Table II, it is easy to see that these expressions are equal, and so 
(9.6) is established for SBC. The proof for DBC is similar. 

VIII. The difference between the natural state S and any completely associated state 
S* is orthogonal to the difference between S and any complementary state S” : 


(S — S*)-(S — 8”) = 0. (9.9) 


Thus the natural state lies on every hypersphere which has for diameter the line joining 
a point of the associated hyperplane to a point of the complementary hyperplane. 

Equation (9.9) is an immediate consequence of (9.1) and (9.6). The statement 
about the hypersphere is verified by writing (9.9) in the form 


[Ss — 3(S* + S")]* = 3(S* — 8”)? (9. 10) 


10. The hypercircle. In (9.1) we located the natural state S on a hypersphere, 
and in (9.5), (9.7) we located it on m+n hyperplanes. This means that the extremity 
of the natural vector $ lies on a hypercircle I’, which is the intersection of the hyper- 
sphere and the hyperplanes. 

The center C of I may be found by starting from the center of the hypersphere 
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(48*), and proceeding through suitable distances in directions normal to the hyper- 
planes until we arrive at a point on all the hyperplanes. Thus, we write 
m n 


C = 38*+ Dad + Dal, (10.1) 


p=1 q=1 


where the coefficients are to be determined from the conditions 


C-I,=0 (p=1,2,---,m), C-I/=S*-I) (¢=1,2,---,n). (10.2) 
By (7.1) and (8.3), we find 
a, = —4S8*-1,, a, = 38*-I;, (10.3) 
and so the center C of the hypercircle T is at 
1 = , ay = La} wh 
C= - |s* — >1,(S*-1,) + DKis+-1) |, (10.4) 
« p=1 q=1 
We note that the magnitude of C is given by 
1 m ; n P 
C? = — & — > (S*-I,)?+3>> (s*-1)*]. (10.5) 
+ p=1 q=1 
We find the radius R of T from 
R? = (S — C)? 
= §-(S — S*) + S-(S* — 2C) + C’. (10.6) 
The first term vanishes by (9.1), and we easily find from (10.4) and (10.5) 
1 m , n ? 
R? = — |S — >) (S*-I,)?- s*x)*]. (10.7) 
4 p=1 q=1 


To sum up: 

IX. The natural state S is located on a hypercircle T with center C as given by (10.4) 
and radius R as given by (10.7). Using X as a current position vector, we may write the 
equations of I in the form 


X-(X—S*)=0, X-I,=0 (p=1,2,---+,m), 


aA v7 (10.8) 
X-I, = S*-I, (¢=1,2,---,%). 


Let us now approach the hypercircle in a different way. 

Points on the associated hyperplane L% satisfy some of the conditions imposed 
on §, and points on the complementary hyperplane L,’ satisfy others. If the two 
hyperplanes had a common point, then (as remarked earlier) that common point 
would be the natural state S. Unless we are particularly lucky in choosing the states 
S*,S/,S/’, this common point will not occur, and then it seems appropriate to seek 
those points, V* on L% and V,/’ on L,/’ such that the distance | VE—V,i’ | is as small 
as possible. We may call these two points the vertices of the hyperplanes. 

The general point on L> is 


Ut ='S* + >> dT}, (10.9) 


p=1 
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and the general point on L,!’ is 


U’ = > &I,. (10. 10) 
q=1 
The square of the distance between these two points is 
m n 2 
(U* — U")? = (s* +>o-> uty). (10. 11) 
p=1 q=1 


It is easily seen that to minimize this, we must choose 
, yon * , i * ” 
b, = — S*-I,, b, = S*-I,. (10.12) 


Substitution in (10.9) and (10.10) gives for the vertices of the associated and comple- 


mentary hyperplanes respectively 


v, =S*— > 1 (S*-1) 


pal 
V.'= DUS*-L). eee 
mer 
We see from (10.4) that 
C = 1(Ve, + Vi’), (10. 14) 
and from (10.7) 
R? = (Vi, — Vi)? (10.15) 


Thus the points V* and V,/’ are the extremities of a diameter of the hypercircle I. 

Any chord of I may be represented by a vector Y; it follows from (10.8) that every 
such chord satisfies Y-I1/’ =0 (¢=1, 2, - - - m), and hence Y-V,” =0. 

We may sum up as follows: 

X. The vertices of the associated and complementary hyperplanes, as given in (10.13), 
are the extremities of a diameter of the hypercircle T, described in 1X. The position vector 
Vi’ of the vertex of the complementary hyperplane is orthogonal to every chord of T 


(Fig. 6). 
The above geometrical statement suggests the inequalities 


Vz SSSVo, (10.16) 
and a formal proof is easy. Any point X on T satisfies (cf. (10.8)) 

(X —Vs)-(X-—V%)=0, X-V,’=S*-Vy (10. 17) 
and if we seek maxima and minima of X? subject to these conditions, it is easy to see 
that we must have 

X =aVn+6V., (10. 18) 


where a and b are undetermined multipliers. When we substitute from this in (10.17), 
the second equation gives b as a linear expression in a, and then the first equation gives 
a quadratic for a. Thus there are just two solutions, and it is easy to verify from 
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(10.17) that they are (a=1, b=0) and (a=0, b=1), ie. X=V" and X=V,.’. More- 
over, by (10.13), 
Vn =S**— > (S*-I})2, 
side (10.19) 


vi" = > (S*-1,)2, 
q=1 


and so, by (10.7), 
Ve —viretR20, Vezve (10. 20) 


Thus V, is the maximum of X and V,/’ the minimum, when X is any vector satisfy- 
ing (10.17). Thus (10.16) is proved. To sum up: 

XI. The vertex of the associated hyperplane is further from the origin than the vertex 
of the complementary hyperplane, and the distance of the natural state from the origin is 
intermediate between them, so that 


vi'= > (S*-L’)? < Sts s*® — > (S*-1,)? = VE. (10.21) 
q=1 p=1 


These inequalities place lower and upper bounds on the strain energy of the natural state S. 

11. Approximations. If § (a set of six stress components expressed as functions 
of position throughout the body) is suggested as an approximation to the natural 
state S, the error, as defined in (4.14), is 


e=|S-S|, (11.1) 


and the squared error is 
e = (S — S)*. (11.2) 
But we cannot in general calculate this error, because we do not know S. 
The hypercircle T comes.to our assistance. The following is obvious: 
XII. If the center C of the hypercircle T as given by (10.4), namely 


1 aos , , = Mt tr 
C= >|s" a >» I,(S*-I,) + p Wis* 1) | (11.3) 
p=1 q=1 


is taken as an approximation to the natural state S, the error € 1s precisely the radius R 
of T, so that by (10.7) 


de? = 4R? = S*? — 9° (S*-1,)?— Do (S*-I,)?. (11.4) 
p=1 q=l 


Since the distance between any two points on I cannot exceed the length of a 
diameter, we have the following result: 

XIII. If any state on the hypercircle T is taken as an approximation to §, the error 
satisfies 
2R, (11.5) 


IIA 


€ 


where R is as in (10.7). The same inequality is satisfied if we take as an approximation 
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either the vertex of the associated hyperplane, 


Vv, =St— D1(S*-1)), (11.6) 
p=1 
or the vertex of the complementary hyperplane 
Vv. = > (s*-l,). (11.7) 
q=1 


We note that since C does not satisfy the equations of compatibility, there 
is no corresponding displacement. Hence, if we are looking for an approximation 
for the displacement throughout the body, V,/’ is the best approximation to take for 
SBC, and V,* the best for DBC. 

We may add that if we take as an approximation 


S=aV.+0Vv, (¢+b=1,a20,520) (11.8) 
then ¢ satisfies 
2RbS¢S2Ra if ba, 
; (11.9) 
2Ra Se 2Rb if acob. 


Finally, we shall show how to use the hypercircle [ to obtain lower and upper 
bounds for the error of any suggested approximation §. It is a question of finding 
the least and greatest distances of § from the hypercircle with the equations (10.8). 
That means that we are to find the minimum and maximum values of 


(X — S)? 


when the vector X is subject to the conditions (10.8). 
It is clear that we must have 


X —S§ = o(2X —S*) + Del, t+ Dal, (11.10) 
p=1 q=1 
where the coefficients are undetermined multipliers. Equivalently, 
(1 — 2a)X = S—aS*+ Doel, t+ Kak: (11.11) 
p=1 q=1 


Hence, by (10.8), 
(1 — 2a)X-I, = c, + S-I, — aS*-I, = 0, 
(1 — 2a)X-I,’ = c,/+ §-I;’ — aS*-I,’ = (1 — 2a)S*-Ty. aati 
If we now substitute in (11.11) the values of c; , c/’ given by these equations, we get 
(1 — 2a)X = V — 2aC (11.13) 
where C is the center of I’, as in (11.3), and V is defined as 


v=s- > UG-Y) _ AIG — $*)-I’]. (11.14) 


q=l 
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We may interpret (11.13) as follows: The points on I at minimax distances from S$ 
lie on the diameter of [ which passes through V, and these points divide the line join- 
ing C to V in the ratios a: —3. We have still to determine a. 

The constant a is found by the first of (10.8). On substitution from (11.13), it is 
found that a satisfies the quadratic equation 


4a°C-(C — S*) + 2a[V-(S* — C) + C (S* —V)] + V-(V—S*) =0. (11.15) 


We may state the result as follows: 
XIV. The square error of any state S lies between 


(X:—S)* and (X.—S)?, 


where X, and X2 are the vectors given by (11.13) when the roots of (11.15) are substituted 
for a. 

Less precise but simpler bounds on the error of § may be established by using, 
instead of the hypercircle I’, the hypersphere with same center C and the same radius 
R. We first find whether § lies outside or inside this hypersphere by computing 


(S — C)? — R?. (11.16) 
If this is positive, § is outside, if negative, inside. If § is outside, we have 
|S—cC|—-Rses|S—C|4+R. (11.17) 
lf S is inside, we have 
R-|S-—C|seSR+|S-Cl. (11.18) 


12. Example. As an illustration of the preceding relations, we consider the torsion 
of a prism of square cross section. We choose the centroid of one end section as the 
origin of coordinates and let the x3 axis coincide with the axis of the prism. To avoid 
carrying along unessential constants, we consider the prism bounded by the planes 
x= +1, x2= +1, x3=0, x3=1 and assume that the shear modulus G and the angle of 
twist ? equal unity. On the end sections of the prism we have 7;=0, and the tangen- 
tial displacements ug (a=1, 2) are given, viz. 


a = 0 on #3 = 0, 
Ua = — €asXs ON 2X3 = 1, 
where €1 = €22=0, €12= — €2, =1. On the lateral faces of the prism 7;=0 (4=1, 2, 3). 
We thus have a case of DBC. 
A suitable completely associated state is defined by 
Ug = — €apXpX3, uz = 0, (12.1) 
or, with the shear modulus G =1, 
Fa3 = — €apXs,; Eas = Es3 = 0. (S*) 
We have 


1 1 I 1 4 8 
S*2 = f f Fazio agd Xd X2 = f f (4 + x2)dxid x2 = (12 ° 2) 
~1 v —1 ae ee 3 
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A homogeneous associated state is defined by 
tH. = 0, tls = o(%1, Xe). (12.3) 


The corresponding vector S’ is 


Eas = ¢.4 Eag = E33 = 0, (S’) 


1 1 
S”? = f f Y,a0,adX1dX2. (12.4) 
-17 1 


As is easily seen, the “warping function” g must assume values of equal magnitude 
but opposite signs in any two points which are symmetrically situated with respect 


and we have 


to any one of the lines x, =0, x2 =0, x1 =X2, X1= —X2. A function of this type is 
Gp = 0) Xs (%1 — 2), (12.5) 
leading to the vector 
, 2p-2 2p—1 2 2 
Ey3 = x1" x2 [(2p + 1) x1 + (2p — 1)aa], 
: 2p—1 2p—2 2 2 , 
Eos = x; % [(2p — 1) x1 + (2p + 1) a9]. (S;) 


We find 
8-4! «te OtM- le 
[2(p + r) — 3] [2(p +r) — 1] [2(p + 7) + 1] [2(6 +7) + 3] 
and 
ry 


“ es (12.7) 
(2p + 1)(2p + 3) 


S*.S/ = 


A complementary state can be derived from a stress function ¥(%, x2) which van- 
ishes along the contour of the cross section. The corresponding vector S” is 


Eas _ Eas 8; Eas = Ess = 0, (S”) 


1 1 : 
5"? = f f VM ol Xd Xo. (12.8) 
~t 1 


and we have 


As is easily seen, the stress function must be symmetric with respect to x; and x2 
and even in both these variables. Such a function which vanishes on the boundary of 


the cross section is 


Yo = (x1'— 1)(x2'— 1), (12.9) 

leading to the vector 
Ey, = 29(x:'— 1)x2" ; 
2-1 (S,) 


2 
E23 = — 29% (x, — 1). 
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We find 
‘7 v7 2 897s?(q + 5 + 1) 


Ss, 8, = ——--—_--_-_--  _-- (12.10) 
(2q + 1)(2s + 1)[2(¢ + s) — 1][2(¢+ 5) + 1] 
and 
7 25g? 
s*-S, = ——_——__- (12.11) 
(2q + 1)? 


We now proceed to determine the vertices Vand V,,’ in accordance with (10.13). 


We have 


I (S*-I/) = Si (S*-S/)/S7?, (12.12) 
where 
2 25.3 
S*.S/ = —-—)> S{? = — (12..13) 
5-7 


by Eqs. (12.7) and (12.6). Thus, 





/ 
li (S*-l!) = — 7 Si, 
and 
Vi = S* —I/(S*-I/) = S* + — Si. (V;*) 
In a similar manner, we find 
‘7 ff, ad 5 vr v7 
V; = J, (S*-I;) 33 Si. (V;) 


To obtain V¥ and Vi’, we must first orthogonalize S/ and S?’ with respect to 
2 ’ 2 
S/ and S/’, respectively. As is easily seen the vector 


S/ = Ss — S/(S/-S/)/s!? (12.14) 


is orthogonal to S/. Using (12.6), we find 


and 
212 


——-: 12.15) 
3°-7-11 ( 


S¢ 2 = Ss 2 (Si -S3 )?/Si 2 


In a similar manner, a vector orthogonal to S/’ can be derived from Sz’. One finds 


4. 
2 3 =f —tik 


7 Si (S2 ) 


9 
14 


= S:’— S$; (S1'S2)/S)" = S:'— 
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and 
a2 72 Malls 2 213. 31 12.16 
Sr = St — Cras: = Sas (12.16) 
We now determine 


IJ (S*-I/) = $4 (S*-S/)/S22. 
With the use of Eqs. (12.6), (12.7) and (12.15), we find 











(S*-l) = — sy. (12.17) 
28.5 
Similarly, 
1,’ (S*-Ii’) ee (12.18) 
25.31 
Thus, by Eqs. (10.13), 
Vi = Vi+ se S/ (V:*) 
28.5 
and 
vim vig 25 ge vs") 
25.31 


In connection with the evaluation of (10.21), we need the values of (S*-I, )? and 
(S*-I7’)?. One finds 


3.7 


/ 
(S*-I/ )? = (S*-S/)?/S/?2 = = > 0.4148 


11 
(S*-IZ)? = ——— > 0.0023, 
38. 52-7 





9 


«. 


i ee 
(S*-1, )? = = * 2.2222, 


9 


—— > 0.0258. 


(S*-I,)? 
iv * Sa 


II 


Now, S? here equals the product of the torque and the angle of twist. Since the 
latter has been taken as unity, S? represents the torque per unit angle of twist or 
the torsional stiffness. Equation (10.21) thus furnishes the following bounds for the 
torsional stiffness :* 

a) for m=n=1: 2.2222 <S? <§$—0.4148 <2.2519. 

b) for m=n=2: 2.2480 <S?<2.2496; 
the exact value, to four decimal places, is S* = 2.2492. 

* C. Weber [Z. angew. Math. Mech. 11, 244-245 (1931) ] used the variational principles of elasticity 
to obtain lower and upper bounds for the torsional stiffness of a square hollow prism. His basic formula 
corresponds essentially to the form which (10.21) assumes in the case of torsion for m=n=1. 
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Next, let us consider the maximum shearing stress, i.e. the stress component £2; 
at the point x; =1, x2.=0. For the various states considered above the value of this 
stress component is given in the following table: 


State | S* | S | { | Si’ |S?’ | S. 


mw 


S w" 
9 








| 
pevignieawe Catia | 2% 
E3(1, 0) | 1 | 1 | 0 | 2 | 4 | — 11/3 | 22-11/(5-7) 
Accordingly, the stress values corresponding to the vertices of the associated and 
complementary hyperplanes are: 


7 25 
1+ —— = — = 1.389 for V¥F, 
2-3? 18 
25 S22 43 431 
— — ——-— = — = 1.347 for V+, 
18 26-5 3% 320 
5 7 
—-2=—=1.250 for V;, and 
2 4 


é. 343 ’ 
4+ a = US 1.383 for Vo. 





The stress values corresponding to 

Ci = (VF + V1) 
and 

C, = 3(V# + V2) 


are therefore 1.319 and 1.365, respectively; the exact value, to three decimal places, 
is 1.352. 

13. Comparison with other work. S. Bergman* has given an exact solution of the 
general boundary value problem for an isotropic elastic body, the surface stress or 
the surface displacement being assigned. In the vector notation of the present paper, 
his solution may be written 


x 
S= i Z,(Z,°S), 

p=1 
where Z, form an infinite orthonormal sequence of states, each satisfying the equa- 
tions of equilibrium and compatibility. These states are determined explicitly. The 
scalar product Z,-S can be computed when either stress or displacement is given all 
over the bounding surface. A point of Bergman’s function space is a set of three func- 
tions (components of displacement), whereas a point of our function space is a set of 
six functions (components of stress), since we relax the equations of compatibility. 
But a more important difference between Bergman’s work and the present theory 
lies in the fact that he was primarily concerned with an exact solution, whereas we are 
interested in comparatively simple approximate solutions in which the error can be 
easily computed. 

It is interesting to discuss, along the lines of the present paper, the results given in 


* Mathematische Annalen, 98, 248-263 (1928). 
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Courant-Hilbert, “Methoden der mathematischen Physik,” vol. 1, J. Springer, Berlin, 
1931, pp. 228-230. We shall refer to this work as CH. We shall show how their re- 
sults can be obtained without recourse to the calculus of variations. Then we shall 
strengthen these results, and express them as a single inequality. Finally, for the 
boundary conditions considered in CH, we shall show that the natural state lies on 
a certain hypersphere. 

The boundary conditions used by CH are different from those of the present paper. 
In their simpler type of boundary condition, they divide the boundary into two parts, 
lr, and Ts, with 7; assigned on T';, and u; assigned on I's. Body forces are included in 
CH, but we shall here omit them. 

Unless either I; or Tl’; disappears, these boundary conditions are neither SBC nor 
DBC. We shall therefore abandon our previous notation, and use the following sym- 
bols: 

S =natural state. 

S’ =a state satisfying the equations of equilibrium and making T/ = 7; on T. 

S” =a state satisfying the equations of compatibility and making uj’ =u,on To. 
Now 


fc — ¢:;)(Ei; — Ei;)dv 2 0, (13.1) 
or 
S? + S’? — 2(e-E’) = 0, (13.2) 
or, in an obvious extension of the notation of (2.10), 
S? + S’? > 2(u-T’): + 2(u-T")2, (13.3) 
where the subscripts 1 and 2 refer to integration over I; and I respectively. Also 
S? = (u-T), + (u-T)>. (13.4) 
Subtracting twice (13.4) from (13.3) and using the fact that T/ = T; on T), we get 
S’2 — §? > 2(u-T’)2 — 2(u-T)o, (13.5) 
or 
4S? — (u-T)2 S 3S’? — (u-T")o. (13.6) 


This is the same as the principle II of CH. By a similar argument, we obtain 
1S? —_ (u- T); < 45”? = (u” - T)1, (13.7) 


which is the same as the principle I of CH. 
If we add (13.6) and (13.7), and use (13.4), we get 


3S’? + 3S"? — (u-T’)2 — (u"-T)1 2 0, (13.8) 


a result which is of no particular interest. 
Now let us strengthen the preceding results, using (3.9), which yields 


S-S’ s SS’, S-S’ s SS’. (13.9) 
The first of these gives, instead of (13.3), the stronger inequality 
SS’ > (u-T’): + (u-T’)2 = (u-T)1 + (u-T")2. (13.10) 
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Subtracting (13.4), we get 
SS’ — S? => (u-T’)2 — (u-T)s, (13.11) 


or, instead of (13.6), 
S$ — (a-T)s 3 SS’ — (0-Ts (13.12) 


Similarly, instead of (13.7), we get 
S? — (u-T); S SS” — (u"-T)). (13.13) 
When we added (13.6) and (13.7) we got the trivial inequality (13.8), which tells 
us nothing about the natural state. But if we add (13.12) and (13.13), and use (13.4), 
we get 


S? S S(S’ + S") — (u-T’)2 — (u"-T)1, (13.14) 


or 
[S — 4(S’ + S”)]? < 3(S’ +S”)? — (u-T’)2 — (u"-T)1. (13.15) 
This single inequality bounds the strain energy of the natural state below and above, 
for all quantities except S are calculable, and we have 
10S + S") —-RSESS34(S'+S") +R, (13.16) 
where R>0O and 
R? = 4(S’ +S”)? — (u-T’)2 — (u"-T). (13.17) 
But we can do better than this. We have 
S? = (u-T)1 + (u-T)s, 
S-S’ = (u-T’), + (u-T’)s, (13.18) 
S-S’ = (u’-T), + (u’-T)e. 


Il 


Since (u- T),;=(u-T’):, (u- T)2=(u"- T)2, we get 


S? — §.(S' + 8") = — (u-T’)2 — (u"-T)1, (13.19) 
or 
[s mS 1(S’ +. S”) |? — R?. (13.20) 


where R?=3(S’+8")?—(u:T’)o—(u"-T);. Thus the natural state lies on a hypersphere 
with center at 4(S’+S") and radius R. 

As remarked in a footnote to CH, their results apply to more general boundary 
conditions. These more general boundary conditions include our SBC and DBC as 
special cases, and may be described as follows: The boundary surface may be divided 
into four regions: 

Tn: JT), Tq assigned. 


Tyo: Tn), Ua) assigned. 
, (13.21) 
T21: 4%), Te) assigned. 


T22: Ucn), Me) assigned. 
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We define 
= natural state. 
S’=state satisfying equations of equilibrium and all boundary conditions on 
T (n) and T «ty. 
S” =state satisfying equations of compatibility and all boundary conditions on 
Ucn) and U4). 
Then 


~ 
i) 
| 


= (U(nyT (ny) + (Uy T ay) 11 + (ny T ny) 12 + (UT) 22 
+ (tenyT (ny)21 + (cy) T ct) 21 + (teenyT (ny) 22 + (UT) 22, (13.22) 
S-S’ = same expression as in (13.22) but with T changed to 7’, (13.23) 
S-S” = same expression as in (13.22) but with « changed to uw”, (13.24) 
Now we have 
(t4¢n)T (ny) 11 = (u(nyT (ny) 115 (ucT (1 = (uT i) 
(t6¢myT (ny)12 = (te(nyT enya» (tec Tey) 21 = (tery Tey) 21 


09 an v7 
(t(nyT (ny)21 = (Ucn) T (ny) 215 (uc) T ce) 12 = (te ceyT (ty) 12 


(t4¢nyT (ny)22 = (tn T (ny) 225 (me) T cey)22 = (eT (ey) 22- 
So we obtain from (13.22), (13.23), (13.24), 
*-8-6 +9) = - (u"- T”)u ad (u-T)2e sae (t6¢ny T (ny) 21 = (ue) T (ey) 12 
= (u(nyT (ny) 12 - (ueyT (ay) 21- (13.26) 


The right hand side is computable. The result may be written 
[S — 3(S’ + 8”) ]* = R?, (13.27) 
where 
R? = 3(S! + S")? — (w"-T")ur — (4 T’)on — (t6¢ny Tomy) 21 
— (uyT (oy) 12 — (HenyT ry) 12 — (Mey T () 22- (13.28) 


We may sum up as follows: 
XIV. With the boundary conditions as described in (13.21), the natural state lies on 
a hypersphere with center at 


C = }(S'+ 8’) (13.29) 


and radius R as given by (13.28), where S’ is any state satisfying the equations of equilib- 
rium and all boundary conditions on stress, and S" is any state satisfying the equations 
of compatibility and all boundary conditions on displacement. 

Added November 18, 1946: The authors are indebted to Professor A. Weinstein 
for showing a connection between the present method and the method given by 
E. Trefftz.* Trefftz dealt with the Dirichlet problem in a plane (to determine a har- 
monic function with prescribed boundary values). The Ritz method gives an upper 
bound to the Dirichlet integral of the solution u, namely, the Dirichlet integral of a 


function v satisfying the boundary conditions. Trefftz supplied a lower bound by 


* Proc. 2nd International Congress for Applied Mechanics, Zurich, 1926, pp. 131-137. 
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taking a linear combination w of harmonic functions, and choosing the coefficients 
so that the Dirichlet integral of the difference (u—w) is a minimum. It should be 
pointed out that, although the torsion problem is reducible to the Dirichlet problem, 
the torsional rigidity is not given by the Dirichlet integral of a harmonic function. 
Hence, in applying his method to the estimation of torsional rigidity, Trefftz had to 
proceed in an indirect way. 

Let us now consider the Dirichlet problem in the notation of the present paper. 
Let a vector in function space correspond to a function with continuous first deriva- 
tives in a plane region R, bounded by a curve B. Let-S* be a completely associated 
vector; this means a function u* satisfying the boundary conditions. Let S” be a com- 
plementary vector; this means a harmonic function u”. Let S be the solution-function u. 
Let the scalar product of two vectors be defined by the Dirichlet integral 


Ou; aes Ou, OU 
S.-S» “SSG —— — }dzdy. 
Ox ie Oy dy 
ou* <i ou* du” Ou” 
S*-.S’ = fy — —-+— — )dxdy -f u* — ds, 
Ox Ox Ty oy RB on 


since u” is harmonic. Since § is itself a completely associated vector, we may write 
S for S* and u for u*. But u* =u on B, and so we have 


S*.S”=S-S"” or S’-(S*—S) =0. (1) 


Then 


But S is also a complementary vector, and so we may write S instead of S”. This gives 
S*-S =S’? or S-(S* —S) =0. (2) 


This last equation is the Ritz equation (4) of Trefftz’ paper. 

The equation (1) is not that of Trefftz. However, it tells us that S*—S is per- 
pendicular to S”. From elementary solid geometry it follows that, if we drop per- 
pendiculars from the extremities of S* and S on the line of S”, these perpendiculars 
have a common foot, say with position vector T”=kS”. We have then 

T’-(S — T”) = 0. (3) 
This is identical with Trefftz’ equation (10), in which the right side is actually zero. 
Obviously, T” is that vector of form kS” which minimizes (S—T”)?; thus T” corre- 
sponds to that harmonic function of form ku” which minimizes the Dirichlet integral 


of u—ku", where wu is the solution. 
To apply the above ideas to the torsion problem, we recall that the torsional 


rigidity (for G=1) is 
Ou Ou 
ff (« Tit <9 on) andy, 
R Ox Oy 


where u is a harmonic function, taking the value 3(x?+y?) on the boundary. If we 
write So for the vector }(x?+-y?), then the torsional rigidity is 


e464 (4) 
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By virtue of (1) and (2) we can build up the theory of the hypercircle, as we did in 
the general elastic problem. We can introduce orthonormal homogeneous associated 


vectors If (b=1, - + - ,m) with u’=0 on B, and orthonormal complementary vectors 
Ij’ (g=1, - ++, m), each harmonic. The center C of the hypercircle is then given by 
(10.4) and its radius by (10.7). We have 
(S—C)I,=0 (p=1,---,m) 
(S) 


(S — C)-I,’=0 (q =1,---,m). 
Thus we may write 
S=C+ RJ, (6) 
where J is a unit vector satisfying 
J-I, =0 (p = 1,--+-+,m) 
JI,=0 (=1,---,n). " 
Thus 
SoS = So-C + RSo-J. (8) 


The maximum and minimum of this expression, as § ranges over the hypercircle, are 


to be found from 


m 


n 
ty! Melt 
So = p Cplp + pa cq1,+<J, (9) 
p=1 q=1 
where the c’s are undetermined scalars. Since I, -I7’ =0, we easily find 


c, = Sol, (p= 1,---,m) 
c¢=Sol, (q=1,--:,n) 


2 es , ~ ” (10) 
¢ =So— DY Gol)? -— DO (So-l)? 
p=1 q=1 
So-J = ¢. 
Hence by (8) 
So-C — RA <= SS S So-C + RA, (11) 


where 
9 ee , = ” aia 
A= s: - 5 (@-r)- Y (ST) 
p=1 q=1 


So we have upper and lower bounds for the torsional rigidity, as given by (4). 
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THE FLOW OF AN IDEAL INCOMPRESSIBLE 
FLUID ABOUT A LENS* 


BY 


MAX SHIFFMAN (New York University) Ann D. C. SPENCER (Stanford University) 


1. Introduction. The irrotational flow produced in an ideal incompressible fluid 
by the motion of an immersed object is described by a velocity potential function 
satisfying certain boundary conditions. For three-dimensional flows the velocity po- 
tential has been determined explicitly only for a very limited category of geometri- 
cally prescribed objects. In the present paper, we shall add to this list by constructing 
the flow about a lenticular object formed by two intersecting spheres, typical examples 
of which are illustrated in Figs. la—c. This lens flow is closely related to the well 
known case of two completely separated spheres, which might from a general point of 
view be considered as two spheres having an imaginary circle of intersection. In fact, 
it is this relationship which forms the point of departure of the present work and 
which will be clarified herein. 





Fic. la Fic. 1b Fic. ic 


We shall consider, then, a lens formed by two intersecting spheres submerged in 
an ideal incompressible fluid and moving at unit velocity in the direction of its axis of 
symmetry. The velocity potential for the analogous two-dimensional problem could 
be easily determined by the use of conformal mapping but this is of no avail for three 
dimensions. Rather, we shall continue the potential function ® analytically as far as 
possible across the two faces of the lens, and determine its singularities. The extended 
function ® is then defined over a multi-sheeted— in general infinite-sheeted—Riemann 
space with the sharp rim of the lens, the circle of intersection of the two spheres, as 
the branch line around which the Riemann space winds. By the nature of the analytic 
consideration, the singularities are located not only at the inverted points but also at 
the points at infinity in the various sheets. It is then necessary to construct funda- 
mental potential functions defined on the multisheeted Riemann space and having 
prescribed singularities at given points. For this construction we shall use the same 
type of procedure as was introduced by A. Sommerfeld in his classical work on the 
theory of diffraction, where a similar question was considered for reflection across 
plane surfaces. 

* Received Sept. 27, 1946. This paper is essentially an excerpt from a report by the authorsentitled 
“The Force of Impact ona Sphere Striking a Water Surface,” AMP Report 42.1R, Feb. 1945, Appl. Math. 
Group at New York University, prepared for the Office of Scientific Research and Development. The 
application to impact will be published separately. 
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2. Procedure for separated spheres. Let ® be the velocity potential (we use the 
convention that velocity = —grad ®). It satisfies the Laplace equation 


$,. + Pry + %,, = 0, 


vanishes and is regular at ©, and satisfies the following boundary condition: on the 
surface of the lens, the negative of the normal derivative of ® is equal to the compo- 
nent of the velocity of the lens along its normal, or 

Ob 


——= cos X. 
On 


To see how to construct the potential function ®, let us first consider the case in 
which the two spheres forming the lens do not intersect and are thus completely 
separated (Fig. 2). 

A classical procedure for finding the flow due to the vertical 

motion of the two spheres is by the method of images (see [1], p. 
131 or [2], p. 420).* The flow produced by a single sphere alone 
moving in an infinite body of fluid is that created by a dipole of 
suitable strength placed at the center of the sphere. Thus, dipoles 
situated at O and O" in the figure would individually give the 
correct flow if the influence of the other sphere could be neglected. 
But the dipole at O’ destroys the correct boundary condition on 
the sphere K, and likewise for the dipole at O and the sphere K’. 
To remedy this, introduce the image of the dipole at O’ relative 
to K (and likewise for O and K’), namely, a dipole of suitable 
strength situated at the inverse point Q; (and at the inverse 
point Q/ ). But again the dipole at Q/ violates the boundary con- 
dition on K, etc. Indefinite continuation of this process leads to a 
convergent expression for the potential function which satisfies 
the boundary conditions on beth spheres. 

We shall apply a generalization of the above process to the 
lenticular case, when the sphéres intersect. The process as it Fic. 2 
stands leads to a serious difficulty, namely that some of the suc- 
cessive images of the dipoles lie directly in the fluid, thus creating singularities 
inside the fluid which violate the requirements of the problem. For examples, in 
the particular case of Fig. 1a, the very first dipoles at the centers of the respective 
spheres already lie in the fluid. The difficulty may be overcome by using multi- 
valued dipoles, the images of which lie in the various replicas (“sheets”) of a Rie- 
mann space and do not insert singularities in that portion of the space in which 
we are interested. However, there is no difficulty in the special cases when the exterior 
angle between the intersecting spheres has the value +/m where m is an integer. This 
is discussed in App. A. It should of course be mentioned that when the exterior angle 
between the intersecting spheres is greater than 7, as in Figs. 1a, b, the velocity of the 
fluid is infinite on the sharp rim of the lens. 

3. Dipoles. We select the toroidal coordinate system (¢, ¥, 6) as most appropriate 
for our domain. These are formed from bipolar coordinates in the plane by rotation 
about an axis, and are defined by 








* Numbers in brackets refer to the bibliography at the end of the paper. 
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i) 
~I 
NR 


x =Prcosg, y=rsing, r+ iz = ai cot 3(y + io) 
*=x*+y", and r, z can be expressed sepa- 


where x, y, 2 are rectangular coordinates, 


rately in terms of o, y by 








sinh o sin y 
r=a- =v z=a 
s—T s—T 
where 
s = cosh a, rT = cosy. 


(See [2], p. 165.) Because of rotational symmetry, the angle 6 will never enter. In the 
r, z-plane, the curve Y=constant is a circular arc connecting the points +a and —a 
on the r-axis, while o=constant is a circle with these two points as mutually inverse 
points. By letting y vary from — © to + without any identification, there results 


an infinite-fold Riemann space with the circle r =a, z =0 as branch line. 








}* 
¢=0 
P 
ss ge Vr+T a -_ 
Py ys-T Po yw=0 
o = Log PP) 
PPo 
Fic. 3 


The inverse distance between a point on the z-axis g=0 at Y= and an arbitrary 
point (¢, y, @) is the absolute value of the quantity 


sin (Yo/2) / $—T 
a 





P(yo) = — ? ee ; 


This represents a pole of strength +1 at the point ¢=0, Y= Wo. We shall consider a 
dipole q(Wo) of strength a| csc3(Wo/2)| /2 situated at the same point, and pointing 


downwards, namely 
: a® tee, : 
q(vo) = — rs csc* (Yo/2) — p(Yo) 
Oz 


\ 


(s — 7)!/2 cot (yYo/2) 
nent y 3.1) 
[s — cos (y — a | 





a {° — 7)'/2 sin (YW — Yo) 
2 “fs — cos (— Yo) |3/2 


The Stokes’ stream function ¢(Wo) corresponding to g(Wo) is (see [2], chapter XV, 


especially page 414): 
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t(Yo) = 4a* csc* (Yo/2)r2p*(Yo) = nnn (3.2) 
2 s 0/4 0 2(s —_ r)'/2[s — cos (y a Yo) ]2/2 Be 





From this it follows that the potential function 
q(Yo) — g(w') 

has vanishing normal derivative on the spherical cap y=c if 
Yo + vo = 2¢ 

For then t(Wo) —t(Wo ) vanishes on Y=c. 


4. Multi-valued dipoles. The quantities g(£), #(£) remain as pote ntial ‘and stream 
functions respectively of a flow if Wo is replaced by a complex number £. Considering 


co 


a, ¥ temporarily as fixed, g(£) and ¢(&) as functions of £ in the complex £-plane are 
multi-valued and have branch points at the zeros of s—cos (¥—§&), that is at 


=yt io + 2kr, an 6, 21,22 




















§-plane 
! ! 
! ! | 
| }|¥ | 
tf | 
| | 
| | | 
by-2t+ic 4) ytic Ly+2neic 
em 
? t yw-ic t 
| | | 
| 
; | 
i |W | 
| | 
, | | 
Fic. 4 


For convenience, cut the £-plane by vertical lines proceeding from the branch points 
to +i asshown in Fig. 4. The functions g(&), ¢(£) are single-valued in this cut -plane 
and take opposite signs on the opposite sides of each cut. We define [s—cos (y—£) }" 
to be positive when € is real. 
The basic multi-value dipole Q(Wo) situated at a (real) point o=0, Y=yYo on the 
z-axis of our infinite sheeted Riemann space is defined as follows 
1 9) 
Oyo) = — 
2ridw & — Yo 
where the path W of integration is the two-branched path indicated in Fig. 4. The 
function Q(Wo) is a potential function in (¢, ) since differentiations can be performed 
under the integral sign, the path W being temporarily fixed. 
For an n-sheeted Riemann space in which values of the angle y differing by 2nz 





dé (4.1) 
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would be identified, the corresponding n-valued dipole Q(Wo) is 


1 £) Se 
OY) = f HO op E—* (4.2) 
2rid w 


2n 2n 


Since cot (E— Wo)/2n is periodic in &£ of period 2mm, and q(£) is periodic of period 27, 
it is easy to see that Q“ (Wo) is periodic in W of period 2nz. 

We now briefly describe the properties of Q, “9), a similar description applying to 
Q(Wo) (see [3], [4] for a complete analysis, the following being but a bare sketch). 


On account of the pole of the integrand in (4.1) at Wo one sees that 
O(Yo) = q(Wo) 


is regular as c—0, Yo (where ¥o#2jr). Thus Q(Wo) behaves like a dipole near 
¢=0, Y=yYo, but is regular near c=0, P=Y~o+2ka7, R= +1, +2,--- since the in- 
tegrand is regular at these points. 

In addition, there are still the poles of g(&) itself at &=2ka, k=0, +1, +2,-- > 
due to the presence of the term cot £/2 (see (3.1)). One easily sees that 


Yo) - — 
2kr — Wo 

as g—0 and y—2kr, k=0, +1, +2,--- (Wo 2jr). Thus the function O(Wo) is not 
simply a multi-valued dipole at Yo but it has certain non-vanishing values at « 
(o=0, Y =2kr) in the various sheets and has no other singularities. 

(A constant value at « may be considered as a pole at ~. If in the definition of 
q(&) the term cot &/2 were replaced by cot Wo/2, the resulting function in (5.1) would 
be zero at the points at ~; this would be a simple dipole. However, a simple dipole 
would lead to great difficulties, and our definition of g(£) is in this case the correct 
generalization to the multi-sheeted case. This illustrates how careful one must be in 
the choice and analysis of the singularities, with special reference to the behavior at 
infinity. ) 

For the special case in which Wo is a multiple of 27, the integrand in (4.1) has a 
double pole. For example, Q(0) has an integrand in (4.1) which has a double pole at 


£=0, and the residue there is 
a siny 0 $—T i/2 sin y 
-— + af —( ————_ —— — = @— = Z. 
2s—T d&\s — cos (WW — &) t=0 S$—T 


Q(0) —-z-0 as «o—0, yO (4.3) 


Therefore 


and so Q(0) behaves like z near the point at ~ in therfirst sheet. A similar analysis 
applies to Q(2kr). 

In the n-sheeted case, it is of interest to consider the special case »=1 where we 
are in the original single-sheeted space. The integrand in (4.2) is then of period 27, 
and we have 


QO (Wo) =s q(Wo) = a/2 cot (Wo/2), Yo x 0, (4.4) 


so Q")(Wo) differs from an ordinary dipole g(Wo) by a constant. If in (4.4) we now let 
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Yo—0, we obtain 
VW) >2z as wo. 
It may also be verified directly that 
QM (0) = z. (4.5) 
5. Inversion of multi-valued dipoles. It is well known that in ordinary space the 


image of a dipole with respect to a sphere is another dipole situated at the inverse 
point. In this section the theorem will be extended to multi-sheeted spaces. 


THEOREM 5.1. The function Q(~o) —Q(Wd ) has vanishing normal derivative on the 
spherécal cap ~ =constant =c if Pot Wd =2c. 

Proof. The stream function T(Wo) corresponding to Q(Wpo) is 

1 t(é) 
T(Wo) = —| eae (5.1) 
2ridw § — Yo 

where ¢(£) is the stream function corresponding to q(£), obtained from (3.2) by replac- 
ing Wo by &. It is required to show merely that T(Wo) — T(¥¢ ) is constant on Y=c. 

Make the substitution £-y=¢. Then 





a a? sinh? ¢ (5.2) 
O= 2 (s — r)'/*2(s — cos £)3/? , 
1 i(£) 
T(Yo) = —f —__—_———d 5.3) 
W 2ridw S+yYP— Wo ; ( 


where W’ is the path in the ¢-plane shown in Fig. 5. By the symmetry of the path W’, 
the lower path W{ is obtained from the upper path W,/ by replacing ¢ by —f¢. The 











| | Wy §-plane | 
| ' | 
| | | 
| | l 
4 U ie 4 
' (~) li ' Whe Why + WH 
| . 
| aty 
198 | 
| im ! 
' | 
Fic. 5 


function ¢(£) is by (5.2) an even function of £, and (5.3) can be written in the form 


1 1 1 
T(Yo) = f t(£) (—. —— + ——— )ar. (5.4) 
2rid wry f+ Ww— yw) &— W— Wo) 
Likewise for 7(W¢ ). Evaluating at Y=c, and noting that the hypothesis pot Wo =2c 
can be written c—Wo = —(c—yW¢ ), eq. (5.4) shows immediately that 7(Wo) —T(We ) =0 
when ~=c, and the theorem is proved. 
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THEOREM 5.2. The function Q(Wo)+Q(—Wo) vanishes on Y=0. In particular Q(0) 
vanishes on y=0. 


Proof. For Y =0 the function g(£) is an odd function of — by (3.1). The same argu- 
ment as above yields Q(Wo) = —Q(—yo) for y=0. 

6. Construction of the potential function for a symmetrical lens. We are now pre- 
pared to apply the method of images described in Sec. 2. For the sake of convenience 
we shall first consider the case of a symmetrical lens, the two spherical portions hav- 
ing the equation Y=8 and Y= —8 where the exterior angle between the two inter- 
secting spheres is 28. The radius of the circle of intersection of the spheres is a, the 
semi-thickness of the lens is B, and the radius of the two symmetrical spheres is A, 
as indicated in Fig. 6. The motion of the lens is in the negative z-direction. 

Let ® be the desired potential function, 


and set 
V=@-—2z. (6.1) 


Then V is a potential satisfying the following 


conditions 





f OV 
4 ——=0 on Yy=8 (6.2) 
/ on 
OV 
——=0 on y=-8 (6.3) 
On 
and 
V behaves like —z at ©. (6.4) 


The potential function V corresponds to the flow about a fixed lens, the fluid having 
unit velocity at o. 

We shall construct V by the method of images, using the infinite sheeted Riemann 
space described by OSa0S ~, —~ <<. (It may be pointed out that if the angle 
8 =n7/2q where n, q are integers, it would be sufficient to use an n-sheeted Riemann 
space in place of the infinite sheeted one.) First start with —Q(0) which has the cor- 
rect behavior (6.4) as a0, y-0 (by (4.3)). To satisfy (6.2), one must add Q(28) 
by theorem (5.1); and to satisfy (6.3) one must add Q(—28). But Q(—28) destroys 
the correct boundary condition on Y=8, and likewise for Q(28) on the spherical cap 
y=-— £6. To remedy this one must add —Q(48) by Theorem 5.1, and likewise 
—Q(—48). Continuing this process indefinitely, one obtains 


V= >> (— 1)*0(2k8). (6.5) 
k=—w 
The convergence of this sum follows directly from its evaluation below. 
We may write (6.5) in the forms 


oo 


V = > (- 1)**[Q(— 248) — Q(2(k + 1)8)] (6.6) 


k=0 
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oo 


> (= 1)**[Q(2k8) — Q(— 2(k + 1)8)]. (6.7) 


k=0 


lI 


By virtue of Theorem 5.1, eqs. (6.6) and (6.7) show that (6.2) and (6.3) are satisfied 
respectively. From the appearance of —Q(0) in (6.5) we see that (6.4) is also satisfied. 
Substituting from (4.1) we obtain 


1 a 1 
v=-—f @® DO (-9 —_ a 
WwW 


2Qri k=—2 §— 2kB 
or 
1 wt 
V=- f o(8) — csc — dé. (6.8) 
2rid yw 28 28 


For the desired potential function ®, we have from (6.1) and (4.5), 


1 T wé 1 é 
o= — f q(é) — CSc - cot dé. (6.9) 
2r1 Ww 2B 28 2 


Formulas (6.8) or (6.9) are the final expressions for V or ® respectively. 
The integrals may be expressed in real form by integrating by parts and then mov- 
ing the path W into the cuts which it surrounds. The results are 


A sin 1p g we T me we 
V= af a/ —— ie — tm [cor csc — + — csc — cot — | ds, 
V/ cosh u—s 2 28 B 28 28 


where £=~+7u, and an analogous expression for ®. 
8 Virtual mass of the fluid. Tle kinetic energy of tke fluid is (17/2) X (velocity 
of object)”, where 


O® 
= eff fi grad ®)*dxdydz = -ff $ — dS, 
lens On 


and p is the density of the fluid, is a normal directed into the fluid. The contribution 
due to a surface integral over a surface which approaches © vanishes. Using the 
boundary condition for —8®@/dn on the surface of the lens, 


M = eff ® cos AdS, EB 
lens 


where \ is the angle between the normal and the direction of motion, indicated in 
Fig. 6. This quantity M is called the “virtual mass” or the “added mass” of the fluid. 
The same formula (7.1) is obtained if the virtual mass is based on the vertically down- 


bo 


ward momentum of the fluid. 
The virtual mass M is more easily calculated by obtaining a different form for 
(7.1).* Formula (7.1) can be written (motion is in the negative z-direction) 


Oz 
u=-off # — dS, 
lens On 


* The longer calculation required in (7.1) was carried out in the original OSRD report with the same 


result. 
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since —0z/0n=cos X. Green’s formula, 


ov Ou 
ff u— —v—IJ]dS = 0, 
on on 


can be applied to the two potential functions ®, z yielding 


O® Oz O® 
M = -ef f z as +f f p — - + )as (7 2) 
ioe ON ms On on 


where the last integral is taken over a surface which approaches «, for example a 
sphere of radius R where Ro. At , ® behaves like a dipole and we may set 

@~D- cos0/R* at (7.3) 
where (cos @)/R? is a unit dipole at th origin pointing in the direction of motion of 
the lens, and D is the strength of the dipole equivalent ta © to ®. Then on a sphere 
of radius R, z= —Rcos 0, and a simple calculation gives 4rpD as the value of the sec- 
ond integral in (7.2). 

For the first integral in (7.2) over the surface of the object, we have 
O® 


— — =cosaA 
On 


OP 
ea z— dS = p z cos AdS = — p-volume of lens 
lens ON 


since cos AdS = +dxdy, the sign depending on the sign of cos X. Substitution of both 
evaluations into (7.2) yields the final formula, valid for any symmetrical object mov- 


and 


ing in an infinite fluid :* 
M = 4rpD — mass of fluid displaced by object (7.4) 
where D is the “dipole strength” at © of the potential function ® as defined in (7.3). 
It remains to evaluate the dipole strength D for the particular lens flow. The 


dipole g(&) from its very construction has strength 
a®| csc? (€/2) | /2 (7.5) 


where | csc? (£/2)| is the piecewise analytic function of the complex variable & equal 
to +csc? (£/2), the sign being determined in such a way that| csc? (&/2) | is positive for 
real & (there are vertical cuts extending from 2n7—10 to 2n7+i). This may also 
be verified directly from (3.1) or (3.2). 

The dipole strength of ® can now be obtained from (6.9). The integrand has no 
singularity near =0 or near £=y (we are considering Y—0), so that the path W may 
be deformed into the path W” indicated in Fig. 7a. Inserting the dipole strength 
(7.5) of g(&), one obtains for the dipole strength D of ®at «: 


1 a’ J ee wt 1 é . 
D=- —f —| esc® (€/2) | csc — — — cot dé (7.6) 
2rid wr 2 28 28 «2 2 


* The same formula applies for any (asymmetrical) object but then D is the component of the dipole 


strength at © in the direction of motion. 
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where now the diagram in the &-plane is indicated by Fig. 7b. The second term in the 
integrand involving cot £/2 can be integrated explicitly and has the value zero. Fur- 
ther, the path W” reduces to Wj’ by replacing & by —£. Making the change of variable 
7r&/2B =z, we have 


a’ Bz 
D =— csc* — csc 2dz (7.7) 
2rid 1 7 
where the path Wj’ is henceforth denoted by Z and is a vertical line slightly to the 
left of the imaginary axis. 
The virtual mass can now be obtained from. (7.4), (7.7) and 


mass of displaced fluid = $2pA#(2 — cos 8) cos‘ 8/2. (7.8) 


The integral in (7.7) is evalu- 




















ated for angles 8 which are rational 

multiples of 7 in Sec. 9. The results §-plane (-stens 
for the virtual mass as a function A | A Y 

of the thickness of the lens are ex- | | 
hibited in dimensionless form in ; 7 | 

Figs. 8 and 9. In Fig. 8 we have veic | 
plotted M/mpA® as a function of y-ic | ¥" wy Ww", 
B/A, where A is the radius of the | 
spheres bounding the lens and B I | 

is the half thickness of the lens. . 

The case of completely separated | | 





spheres, B/A>2, is included in Fic. 7a Fic. 7b 

the graph (see App. C). It is also 

of interest to emphasize the presentation surface of the lens, and for this purpose we 
have plotted //mpa‘ as a function of B/a in Fig. 9, where a is the radius of the sharp 
circular rim of the lens. For B/a>1, the lens has a figure eight cross section and the 
presentation radius is preferably A instead of a, so that Fig. 8 should be used for this 


range. 


Virtual Mass TrpAjm 
of a Symmetrical Lens 





Fic. 8 
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8. Potential function and virtual mass 
for an arbitrary lens. We construct the 
velocity potential of the flow around an 
asymmetrical lens bounded by the spherical 
caps Y=a—6 and y= —8, where 6B and a—f 
are both positive. We suppose that the lens 
is immersed in an infinite fluid and that it is 
moving downward (i.e. in the negative z-di- 
rection) with unit velocity. The symmetrical 
lens is obtained by taking a= 28, and this is 
the case which has been discussed in the 
main body of the paper. 

The potential function describing the 
flow satisfies conditions (6.1)—(6.4) as in 
& Sec 6, but. in place of (6.2), we have 


Virtual Mass Coefficient K 
M = Trpa?-K 





——=0 on Y=a-—-68 andon y= —-8. (8.1) 
on 
In order to illustrate a remark made in Sec. 6, namely that the m-fold Riemann space 
may be used if the intersection angle, here equal to a, is rational, we suppose that 
a=n7/m. The formulas which we obtain are then valid for all a by continuity. 
The solution is the potential function 


>X [Q™(2ka — 28) — O(— 2ka)] 
sia (8.2) 


YE (2(k + 1)a — 28)— O(— 2ka) | 


V 


ll 
> 
3 


~ 
- 


making use of the fact that Q)(Wo) is periodic in Wo of period 2m7. The first form of 
(8.2) shows by Theorem 5.1 that V has vanishing normal derivative on y = —8, while 
the second form yields the vanishing of the normal 
derivative on Y=a—f. Furthermore, at in the 
first sheet (o=0, Y=0), all the terms in the 
summation are regular except —Q(—2ma) = 
—Q™(—2n7r) = —Q™(0) which behaves at ~ like 
—z+0(1) (compare (4.3)). Substituting from (4.2) 
into (8.2) and evaluating the sums, we obtain the 
formula 


1 T m(é + 28) r 
V= —f g(é) - (cot \ - “— cot = at (8.3) 
W 


21 2a\ 2a 2a 





The desired potential function = V +z. Fic. 10 

Formula (8.3) reduces to the old formula (6.8) in 
the special case when a=28. It is interesting to note that (8.3) is built out of two 
terms of the type 
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ew 1 7 m(E—yo) 
Oyo) = =i q(é) 2a ct — — dé. (8.4) 


According to Theorem 5.1, the image of —Q(0) relative to Y= —8 is Q(—28), while 
its image relative to Y=a—8 is 0(2a—28). But 0(—28) =0(2a—28) since O(Wo) is 


periodic in Wo of period 2a. Thus the potential function 
V = — Q(0) + Q(— 28) (8.5) 


satisfies the boundary conditions on both y= —6 and Y=a—, and has the required 
behavior at «©. It is therefore the solution to our problem. This is exactly (8.3). 

The expression Q(o) in (8.4) might be interpreted as a fundamental dipole in a 
space obtained by identifying angles y which differ by 2a. This would be analogous to 
the two-dimensional case in which V could be extended across y = —8 into the region 
between Y= —6 and Y= —a-—f by inversion and then noting that this extension is 
likewise the extension across Y=a—® if one identifies Y=a—8 with y= —a—B. But 
the simple extension by inversion does not apply in the three-dimensional case but 
rather a more complicated continuation. In fact Q(Wo) in (8.4) is not single valued in 
the space obtained by identifying angles y which differ by 2a. Thus the analogy with 
the two-dimensional case given by formulas (8.4), (8.5) is not exact, but is neverthe- 
less noteworthy. 

The calculation of the virtual mass is similar to that of Sec. 7. The result is 


M = M,— M, (8.6) 


where 


_, mB 1 _ as 7B 7B e 
M, = 2zpa’* sin - | csc®— csc 24csc [| zs + — } — csc{ s — — }pdz (8.7) 
a 2lriJd 1, T a a 


and where 


I B B 
M, = — roa? (2 — cos B) cot — csc? — 
6 2 ? 





+ (2 — cos (a — B)) cot <= csc? va - | (8.8) 


is the mass of the fluid displaced by the lens. 

9. The symmetrical lens. Virtual mass in special cases. (i) For the determination 
of the virtual mass there remains but to evaluate the integral in (7.7), and this can 
be done in terms of residues when the angle 8 is a rational multiple of 7. Suppose that 


-=— (9.1) 


where » and gare relatively prime integers and n <q. The integral to be evaluated is 
then 


1 nz 
J =- —f csc* — csc 2dz (9.2) 
2r1 L q 


where we repeat that L is a vertical line slightly to the left of the imaginary axis. 
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Let 6 be a fixed positive number and introduce the factor e~* in the integrand in 


(9.2). Then 


1 nz 
—f csc* — csc ze~**dz = — sum of residues to the right of L 
2r1 L q 
(9.3) 
= a 1 sum of residues 
1 — (— 1)"t+ae-*a © between 0 and gz. 
Calculating the residues between 0 and gz, letting 6-0 and observing that the coeffi- 
cient of 1/6 must vanish, we obtain two different results according as n+q is odd or 


even. If n+q is odd, 


: i ee nar 
I - = 1)? 1 eset 2 


2 j=1 q 
gq *! gq? d? ; 
d~(- yt fese s+ - (csc | (9.4) 
2n k=1 n* dz? z=kqr/n 


and if m+ q is even, 


1 . q — jnr 
[= “| (= + ) + >> (— 1)*j csc? — 
qrL\6n* 2n jul q 
q n—1 f q’ d? \ 
+ — (— 1)*<z2cscz + - (g csc 2) : (9.5) 
2n X ] n? dz? — f, iaiede 


The computations based on these formulas were used in plotting the graph of Fig. 8. 
It is interesting to note that formulas (9.4), (9.5) involve angles which are multi- 
ples of gr/n as well as those which are multiples of naz/g. These trigonometric sums, 
which would normally depend in a highly discontinuous manner on the angle 8 = n7/q, 
do indeed depend continuously since they are special cases of the integral (9.2). 
For the purpose of checking, it is worth while to calculate the virtual mass in 
certain easy special cases where it can be calculated by other means. 
(ii) If B=7/2, eq. (9.4) for n=1, g=2 gives J=1/2, or from (7.4), (7.7), (7.8) 


M = 2xpA’. 


This trivial case corresponds to the flow about a sphere in an infinite fluid, and the 


virtual mass is then known to be as stated. 
(iii) If we let B—0, we obtain from (7.7) and the fact that a=A sin 8B, 


2) 


3A’ CSC 2 
iA -f —— dz = A* > (— 11/78, 
L . j=l 
Substituting in (7.4), (7.7), (7.8), the virtual mass of the fluid, when the two spheres 


are tangent (G=0), is 


M = tnpa*| Be (— 11/72 - 2/3. (9.6) 


This value agrees with the limit obtained when completely separated spheres ap- 
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proach tangency (see App. C) and so is a consequence of known formulas (see [2], 
p. 447). 

(iv) Taking B=7, the virtual mass M is obviously zero. However if we keep a 
(instead of A) as fixed, then we obtain the virtual mass for a circular disc of radius a. 
Taking B = in (7.4), (7.7), (7.8), we obtain 


4rpa* 
M = - csct 2dz = $pa’*. 
2nt L 


This is also a known result (see [2], p. 428). 
10. The symmetrical lens. Derivative of the virtual mass. For the application of 
the lens flow to the impact of a sphere on a fluid surface, it is necessary to calculate 





dm 

— where m= M/nxpA*, b= B/A. 

db 
This can be done by differentiating formulas (7.4), (7.7), (7.8), leading to another 
trigonometric contour integral. This integral can then be evaluated for B=nm/q as 
in the preceding section, with two cases appearing, according as m+q is odd or even. 
We shall not write down these formulas since they are complicated and are similar 
to those in the preceding section. Just the remark should be made that in equating the 
coefficient of 1/6 to zero one obtains evaluations of certain trigonometric sums, in- 
volving angles which are multiples of gr/n and nz/g, which appear at first glance 
quite startling. 

On the basis of these formulas, a graph of dm/db was prepared and is shown in 

Fig. 11. 





on 
db 
- Derivative dm 
: of the 
Virtual Mass Tp AJm 
6 
4 
on 
0 4 8 1.2 ©6116 2.0 2.4 2.8 b 


Fie. 11 


APPENDIX A 


Figure eight lens with exterior angle a=2/m. The elementary process of succes- 
sive images of dipoles fails in general without the introduction of multi-sheeted spaces. 
But this method succeeds for the special case of a lens with exterior angle a=2/m, 
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in which case the procedure described in Sec. 2 terminates after m steps. On the mth 
step, the same dipole corrects the boundary conditions on both spheres K and K’ 
simultaneously, and the total number of dipoles is 2m —1. An example is a figure eight 
lens with a=7/2, which has three dipoles situated at the 
points O, O’, Q indicated in Fig. 12. (By inversion relative 
to K or to K’, the circle OR O’S inverts into the straight 
line ROS. Thus Q is the inverse of O’ relative to K and also 
of O relative to K’.) The dipoles are directed vertically 
downward and their strengths at O, O’ and Q are 

A/2, A’8/2 and — A%A’3/2(00’)§ = — a3/2 respectively. 


2, 2 

For the general case in which a=z/m, it is convenient 
to use the angular coordinate y. Let the equations of the 
two spherical caps K and K’ be y=—a/2+y and 
YW =a/2+v7 respectively, where —a/2<y<a/2. The dipoles 





are situated on the axis of symmetry at the points Fic. 12 
P,, w=—-(w+1at [i+(—-1)’]y, v=0,---,m—-1 
Pi; w= (t+Dat [1 + (-— 1)’}y, v=0,---,m—1. 


The points P,,. and P,- are identical, since Ya1—Ym1=22. The strengths of the 
dipoles at P, and P; are (—1)’ \SP?/2, where S is a point on the sharp circular rim 
of intersection of the two spherical caps and SP, is the distance from P, to S. 

The expression for the velocity potential is now easily written down in terms of 
these (2m—1) dipoles, and the virtual mass obtained from formula (7.4). 

APPENDIX B 

Spherical shell and figure eight lens with exterior angle a=27/m. In the special 
case when the exterior angle a = 21/m, the velocity potential V (or ®) can be expressed 
in terms of elementary functions. According to Sec. 8 all that is required is a two- 
sheeted Riemann space and the fundamental dipoles can then be explicitly integrated. 

It is much simpler to obtain an expression for the Stokes’ stream function instead 
of the velocity potential. Let 7 (Wo), S, x be the Stokes’ stream functions conjugate 
to O™ (Yo), V, ® respectively. Corresponding to (4.2), (8.3) and ®= V+z we have 


1 t(é) .~? 
T™ (Yo) = — f — cot —— dé (B.1) 
2rid w 2n 2n 
1 mt(£) m(é + 28) mi 
S= —| —~ | eo wl cot | dé (B. 2) 
2ridw 2a 2a 2a 
x=S+7°/2 (B. 3) 


where ¢(£) is obtained from ¢(Wo) in (3.2) by replacing Wo by &, and r?/2 is the stream 
function conjugate to z. In any numerical work these have a much simpler form than 
the corresponding formulas for the potential. 

In the integrand the substitution £=y~+¢ can be made, and the lower half W/ 
of the path of integration can be reduced to the upper half W,/ by the substitution 


¢——f¢. There results 
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(B.4) 





Tin( a’ sinh? o f 1 sin (¢/n)dt 
(n) 5) = — . anon ° 
v w’, (s — cos £)8/2 cos [(W — yo)/n] — cos (¢/n) 


4rni(s — r)*/? 
Consider now the case n= 2. Make the substitution cos (¢/2) =cosh (4/2) cosh 7 
and evaluate the integral. The result is 


| \ (2/2) + cos [(¥ — _ "T 





2 
T® 0 = aa a —_—_—$<——$. $$ ne 
oe = he T aia cosh (¢/2) — cos [(W — Wo)/2] 


21/2q? sinh? (¢/2) cos [(y rag ¥o)/2] (B.5) 


m(s — r)'/2[s — cos (W — yo) | 





where {(Wo) is the fundamental dipole in (3.2). 

The arc tangent factor in (B.5) illustrates the difference in behavior of T® (Wo) 
between Yo, 0 and y-yo+27, c—0. In the former situation, the factor ap- 
proaches 1 while in the latter it approaches 0. The second term in (B.5) just cancels 
the remaining singularity, with the result that T@(W»o) —t(Wo) is regular as Yo, 
g—0 while 7 (Wo) itself is regular as PN-yo+27, c—0. 

We may now obtain the stream function for the flow about a figure eight lens with 
a=2n/m. For the stream function S conjugate to V we have a formula analogous to 
(8.2): 


S = >} [T(4kx/m — 28) — T(— 4kr/m)]. (B.6) 
k=1 


The case m=1 is especially noteworthy since this corresponds to the flow about 
spherical shell (exterior angle 27): 


S = T®(— 28) — T@ (0), a = 2m (spherical shell). (B.7) 
For the stream function x, conjugate to ®, we have 

x = T?(— 28) + T@ (2x),  a@ = 2m (spherical shell) (B.8) 
since 72/2 =T® (0) =T@(0)+T7 (27). The flow about a spherical shell has also been 


obtained by Basset [5] using spherical coordinates. 
For the sake of completeness, we include here the formula for Q® (Wo). It is 


a (s—r)'/2 sin [(W — po) /2] 
Yo) = g(Yo)f 0) + — ——_— ee 
VOW (Yo) f(y 7/2 =~ anes <a 


— = esc (Yo/2)f(0) + — tan (o/4), 


where f(Wo) is the arc tangent factor of ¢(Wo) in (B.5), 
a P 
f(¥o) =—arctan[ |]. 
7 


The potential function arising from a pole (instead of a dipole) of suitable strength 


has the form 
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eI. .koh 0k 
[s — cos (W — yo) }!/2 it 


In particular if the pole is at infinity (Yo=0), its form is 
f(0). 


This last expression is a potential function, equal to 1 at infinity in the first sheet 
(y =0), equal to 0 at infinity in the second sheet (¥ =27) and having no singularities 
in the finite. 

APPENDIX C 

Completely separated spheres. Although the case of completely separated spheres 
can be treated by the method of images described in Section 2, and the virtual mass 
computed from (7.4), nevertheless it is of interest to trace a closer relationship to the 
case of intersecting spheres. It will be shown that the integral formulas for lens flow 
carry over to separated spheres by merely allowing the coordinates o, Y and the geo- 
metrical quantities a, B, a to be pure imaginary. 

For simplicity we shall discuss only the symmetrical case, the extension to the 
asymmetrical case being clear. Let the two spheres be of radius A, as indicated in Fig. 
13, moving at unit velocity in the direction of their line of centers. The distance be- 
tween the centers of the spheres is 2(B—A). As in the lens case, set 

B = A(1 + cos 8), a=Asing tC) 

where now @ and a are pure imaginary. The points 

z marked a and —a in the diagram are mutually 
inverse with respect to both spheres. As in the lens 

case, the dipoles describing the flow are situated 


rt at the points r=0,z=a cot kB fork=+1, +2,---> 
with strengths (— 1)#+1 a*® csc* (kB) /2| respec- 
tively.* 


The virtual mass M can be obtained from 
r formula (8.4) with the result 


M = sin 3 B 


A 
< 8 
= m = —— = 4)) (— 1)*#1* ——_-—, (C.2) 
A apA* ram sin® kB 3 
B . —— ° 
; where it should be recalled that 6 is imaginary 
| and so the series is rapidly convergent. In par- 
t ticular, if the spheres approach tangency, B=0 


and the result (9.6) is obtained. 
As B varies from 7 to 0 along the real axis, the 
Fic. 13 spheres first form a thin lens which thickens until 
the spheres are tangent. Then, as 6 varies from 0 
to ic along the imaginary axis, the spheres separate and move farther and farther 








* We remark that in any problem involving successive inversions relative to two equal spheres, the 
introduction of 8 by means of (C.1) is a simplifying device. A similar device applies to any two spheres, 
and will simplify such formulas as in [1], pp. 131-133 and [2], p. 420. This remark also applies to the elec- 
trostatic analogue, as in [3], pp. 642-649. 
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apart. The virtual mass m and its derivative dm/dg vary continuously, but we see 
from Fig. 11 that d?m/db* has a discontinuity at B=0 (6b=1-++cos B=2). 

We may also obtain the potential function in a form analogous to (6.8) or (6.9). 
Allow the coordinates o, y as well as a, B to be pure imaginary. (More generally, com- 
plex values of o, y, a, 8 may also be selected.) Let r, z be connected with a, y by the 
same relations, 


a sinh o asin yp 








' cosh o — cos ry "cosh ¢ — cos y 

The coordinates r, z are still real when a, y, a are all pure imaginary. In order to have 
r=0, we suppose that 0S Imam while Imy varies from — © to +. The locus 
y= —8 is the sphere of radius A and center at r=0, z= —A cos B. The coordinates 
Imo, Imw and the angle @ of rotation form the so-called bipolar coordinate system 
in space. 

The function defined by (3.1) remains a potential function even though oa, y are 
pure imaginary. The branch points of g(&) as a function of the complex variable & 
are still located at ~—io+2kr, ~+io+2kz7 but all of these now lie on the line 
Im £=Imy. The &-plane can be cut along the straight line segments connecting the 
two adjacent branch points y—io+2km and y+io0+2(k+1)z (see Fig. 14). Select the 
path W as indicated in Fig. 14; W can also be considered as surrounding all the indi- 


vidual cuts as well. 
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The potential function V (or ®) can now be constructed as in the lens case, with 


the same form as in (6.8): 


y= @ <ce Sa (C.3) 
= gé 2B csc 2B 2 


2nrid #F 


The formulas for the virtual mass M is the same as (7.4), (7.7) except that the mass 
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of the fluid displaced is merely §1pA*. The integral in (7.7) can then be evaluated by 
residues for 8 pure imaginary, with the result already obtained in (C.2). 
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VIBRATIONAL PROBLEMS IN ELLIPTICAL COORDINATES* 


BY 
N. W. McLACHLAN 
(University of London) 


1. Introduction. We commence with the differential equation 
d°*Fo d°Fo df 
—, a ie K Sipeaie = 


Ox? 0 y” Ot? 


0, (1.1) 


K being a constant. Writing ()=e‘¢(x, y) in (1.1), it becomes 
ay ag 
—+— + & = 0, (1.2) 
dx2 dy? 
with kj =w*K. This is the familiar two-dimensional wave equation for sinusoidal 
time-displacement, expressed in rectangular coordinates. The ordinary differential 
equations used herein, into which (1.2), expressed in elliptical coordinates, may be 
decomposed, are identical in form with those obtained by Mathieu! when he solved 
the problem of the elliptical membrane. Their preferred canonical forms are 


d*y 
- + (a — 2q cos 2z)y = 0, (1.3) 
a2* 

and 
d*y 
rt (a — 2q cosh 2z)y = 0, (1.4) 
dz? 


with g>0. In dealing with an elliptical plate, these equations are used with q20. 
(1.3) is known as the ordinary Mathieu equation, and (1.4) as the modified Mathieu 
equation. They are derivable from each other by the substitution +7 for z. Recently 
lists of solutions of (1.3), (1.4), and the corresponding equations for g<0, have been 
published.?* By aid of these we shall derive formal solutions pertaining to the vi- 
brational modes of: 

(a) A uniform, homogeneous, loss-free stretched membrane in the form of an 
elliptical ring ;5 

(b) Water in a lake of uniform depth whose plan view is an elliptical ring; 

(c) A uniform, homogeneous, loss-free, elastic elliptical plate; 





* Received Oct. 28, 1946. 
1E, Mathieu, Mémoire sur le mouvement vibratoire d'une membrane de forme elliptique, Jour. de Math. 
Pures et Appliquées. 13, 137 (1868). 
2W. G. Bickley, and N. W. McLachlan, Mathieu functions of integral order and their tabulation, 
MTAC, 2, 1, 1946. 
3N.W. McLachlan, Mathieu functions and their classification, Jour. Math. Phys., 25, 209, 1946; 
Theory and application of Mathieu functions, Oxford Press, New York, 1947. 
4 J. A. Stratton, P. M. Morse, L. J. Chu, and A. Hutner, Elliptic cylinder and spheroidal wave func- 


tions, Wiley, & Sons New York, 1942. 
5 An elliptical ring is that portion of the ellipse enclosed by the inner and outer confocal bounding 


ellipses. 
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(d) A uniform, homogeneous, loss-free, elastic elliptical ring plate. 

In all cases, as the outer bounding ellipse tends to a circle, the formulae degenerate 
to those already known for a circular boundary. 

2. Elliptical ring membrane. In this case kj =w*pi/r, where p;=mass of membrane 
per unit area, 7 =radial tension per unit arc length,® both of the latter being constant. 
To derive the appropriate formal solutions, we first transform (1.1) to elliptical 
coordinates,’ where 


x = hcosh é cos n, y = hsinh é sin n, (2.1) 


2h being the interfocal distance. Then if ¢ is the displacement of a point (&, 7) on 
membrane, we get 
4 07¢ 
+ + 3k*(cosh 2 — cos 2n)f = 0. (2. 
Og On? 
with k?=1kih? =wh’p,/4r. 
3. Solutions of (2.2). The appropriate type take the form 


o(E, n) = x(&)W(n) (3.1) 


- 


where W and x are functions of , — alone, respectively, which satisfy the ordinary 


equations 
d*y 
- + (a@ — 2q cos 2n)y = 0, (3.2) 
dn? 

and 
d*x 
; — — (a — 2q cosh 2§)x = 0. (3.3) 
dé 


a is the arbitrary separation constant, while g=k?, k>0. These equations have the 
same forms as (1.3), (1.4), and we designate their first and second linearly independent 


solutions by 
¥iln,q), vel, 9); (3.4) 
xilé, q), x2(é, q). (3.5) 


As shown in the papers mentioned in Footnotes 2, 3, 4 these solutions have a plu- 
rality of forms, so the next step is to select those appropriate to the problem of the 
ring membrane. 

4. Physical conditions. Consider the maximum displacement of all points on a 
confocal ellipse on the membrane. If we start at 7=0 and move round the ellipse 
counterclockwise, the maximum displacement varies continuously. If may be re- 
peated between »=7 and »=2r, but it is always repeated at interval 27. Thus ¢ is 
single-valued and periodic in the coordinate n, the period being either a or 27, ac- 


®° To obtain uniform tension, the membrane may be stretched with uniform radial tension over a 
circular frame, the elliptical rings then being clamped in position within the frame. The plane of the 
membrane, when at rest, is horizontal. 

7 A simple method of transformation is exemplified in Phil. Mag. 36, 600, (1945). 
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cording to the mode of vibration. Consequently y(n) and W2(n) must be functions 


such that either 
¥(n, 9) = ¥(n + =, 9), (4.1) 


or 


V(n, g) = W(n + 2r, q). (4.2) 


According to the references in footnotes 2, 3, the only functions which satisfy these 


conditions are, respectively, 


¥i(n, q) = Ceon(n, 9), S€2n+2(n, Q), [aon, bons2],® (4.3) 
and 

Vilm, 9) = Ceonsi(m g), — Se2nza(ms qs [danas bensa],® (4.4) 
n=, 1, 2,--- The functions in (4.3), (4.4) are first solutions of (3.2). Linearly 


independent second solutions are non-periodic and, therefore, inadmissible here. 
There is no physical reason for discriminating between independent solutions of (3.3) 


so we may include both. Then,?# 


Cen(é, q), Sem(é, q), [am, bm], (4.5) 


xi(é, q) 
and 


Feyn(é, q); Geyn(&, q), [@m, bn, (4.6) 


xo(é, q) 
where m=2n, 2n+1, 2n+2, as the case may be. The second alternative solutions 
Fe,,(E, g), Gem(~, g) may be used for x2(&, g), but for reasons stated in the references 
mentioned in Footnote 3, the solutions (4.6) are preferable. 

5. The formal solution. This is to be constructed from (3.1), and (4.3)—(4.6). We 
arrange the xy in groups, the functions in each group having the same characteristic 
number. Thus, introducing the real part of the time factor from Sec. 1 with an arbi- 
trary phase angle, we have the component solutions of order m, corresponding to @m 


and b,, namely, 


om(E, n, 2) [CmCem(E, 9) + FinFevm(E, 9) |cem(n, 9) COS (mt + €m), (5.1) 


and 


[SmS€m(E, 7) + GmGevm(E, 7) |sem(n, @) COS (mt + Em). (5.2) 


Em(é, 0, 2) 


In (5.1), (5.2), Cn, Fm, Sm, Gm, are arbitrary constants determinable from the dis- 
placement and velocity over the surface of the membranal ring at ¢=0, while wm, Gm 
are the respective pulsatances of the mth free modes of vibration, and €,,, é, their 
relative phase angles. w, and @, are different as will be evident from the next section. 
For Ce, ce, Fey, m=0, 1, 2, - - - , while for Se, se, Gey, m=1, 2, - - - . The complete 
solution, without the time factor, maybe written 


¢(t, n) = dX talé ») + DX Enlf, 0). (5.3) 
m=1 


m=0 


8 These are the characteristic numbers corresponding to the Mathieu functions ce, se of integral order 


2n, 2n+1, 2n+2. 
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6. Vibrational modes of ring membrane. We use the boundary conditions 


£=£, ¢=0 at the outer clamp, £=£,, ¢=0 at the inner clamp. For modes of order m 
S$—S0» § Ss I 
corresponding to characteristic number dm, by (5.1) at the outer clamp 


[CmCem(Eo, 9) + FF eym(Eo, 9) |cem(n, g) = 0. (6.1) 


At the inner clamp 


[CmCem(€1, 9) + FmFevm(E1, 9) |cem(n, g) = 0. (6.2) 


These equations are independent of 7 at the clamps, so for C,,, F, to be non-zero, 


we must have 
Cem(Eo, Q)Feym(é1, 9) — Cem(E1, g)Feym(Eo, 9) = 0, (m = 0) (6.3) 


which is the pulsatance equation. Moreover, corresponding to each value of gm.p, 
p=1, 2, 3, - - - for which (6.3) is satisfied, the elliptical ring has a vibrational mode 
of order m and rank p. Since gm,p=i,,p =Wm,ph*pi/47, the pulsatances of the modes 


are given by 


2 
< 


Wm,p = 47Qm.p/ h*pr, (6.4) 
4rqm,p/a7e"pi, (6.5) 


where h =ae, a being the semi-major axis, and e the eccentricity of the outer bounding 
ellipse. 9m,» are the parametric zeros of (6.3), and determine the positions of the nodal 
ellipses, of which there are (p—1) within the clamping rings. The nodal hyperbolae 
are determined for any g=Qm,p» by 

Cém(n, g) = 0, (m > 0). (6.6) 


For the set of modes corresponding to characteristic number ),,, the pulsatance 
equation defining the nodal ellipses is 


Sem(Eo, Q)GeYm(E1, 9) — Sem(E1, G)Geym(Eo, g) = 0, (m > 0), (6.7) 


where ¢=Qm,p. The nodal hyperbolae, for any g=Qm,p, are given by 
Sém(n, g) = 0, (m > 0). (6.8) 
7. Elliptical ring lake. The formal analysis is identical with that for the ring 
membrane, but the condition at the inner and outer boundaries is zero veiocity 


normal to them. Thus if ¢ represents the tide height or vertical displacement of the 
water from its equilibrium level, we have® 


o¢/dE = 0 at E=, andalsoat §& = &. (7.1) 


The conditional equations have the forms (6.1), (6.2), but the functions are replaced 
by their first derivatives with respect to £. Hence by (6.3) the pulsatance equation is 


Cem(Eo, Q)Fevn(é1, 9) — Cem(Es, Q)Fevm(to, 9) = 0, (m 2 0), (7.2) 
q having those values g),,p, p=1, 2, - - - for which (7.2) is satisfied. These determine 


the nodal ellipses, while the nodal hyperbolae for the same values of g are given by 


9H. Jeffreys, Free oscillations of water in an elliptical lake, Proc. Lond. Math. Soc. 23, 455 (1924). 
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(6.6). In the present problem k? = q=w°h?/gd, g being the acceleration due to gravity, 


and d the uniform depth of the lake. 
For the set of modes corresponding to characteristic number b,, we have 


Sem £0, g)Geym(E1, q) apa Sen(é1, q)Geyn(Eo, q) aad 0, (m > 0), (7.3) 


and 
Sen(n, g) = 0, (m > 0), (7.4) 


with g=G,,.p 

8. Elliptical plate. Let p be the density, ¢ the uniform thickness, ¢<1 Poisson’s 
ratio, E the modulus of elasticity, and c?=Et?/12p(1—¢), whose dimensions are 
(velocity).2 Then if ¢ varies as e**‘ in the various vibrational modes, kj =w*/c*, and 
rotation effects are ignored, the differential equation of motion expressed in rectangu- 
lar coordinates, may be shown to be 














0*¢ 04 20%¢ 
a ae 8 _— 
or 
([+5+H\(= +<-#)r=0. (8.2) 
Ox® = §=dy? Ox? = dy? 
If a function ¢; satisfies 
“ — + kits = 0, . (8.3) 
and another function ¢- satisfies 
SS Ww iene (8.4) 


Ox? Oy? 


their sum, with the appropriate arbitrary constants, is a solution of (8.1), (8.2). To 
obtain this solution, (8.3), (8.4) are expressed in elliptical coordinates, thereby 
yielding two equations of the form (2.2), with +? in one and —&? in the other. 
These may be split up into pairs of equations like (3.2), (3.3), one pair being asso- 
ciated with +q and the other with —g. Then omitting the time factor, the solutions 


of order m corresponding to (8.3) are 
(m) 
1 (é, n) — [CmCem(€, q) + F,,Feym(é, q) |cem(n, q); (a _ Gm), (8.5) 
~(m) " hn 
Fi (é 2) = [SmSem(é, g) + GmGeym(E, q) ]sem(n,q), (a = bm). (8.6) 
The solutions corresponding to (8.4) are 
(m) apf ee Py ~ 
fe (é, n) - [CnC em(E, — q) + F,,Fekm(&, _— q) ]cem(n, = q) (8. 7) 
for @=d2,, M=2n; @=Deny;, M=2n+1, and 


3 (m c Apr 
62 ¢¢, n) - [S mSem(E, ns q) + GinGekm(&, = q) |sem(n, sat q) (8.8) 
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for @=bo,, m=2n;a=2n+1, m=2n-+1. In these cases g <0 in (2.2) and the functions 
Fekm(~, —q), Gekm(€, —q) have been used in preference to Feyn(&, —q), Geym(&, — 4), 
since the former are real, while the latter are complex if & is real.* 
On the interfocal line of the ellipse £=0. In crossing this line orthogonally from 
(0, 7) to (0, —7), we must have 
(a) continuity of displacement, so ¢(0, 7) =¢(0, —7n), and 
(b) continuity of gradient, so 
— [r(6, n) Jeo = — — [t&, — Deo 
0& re) 


ctr 


Remembering that ¢(&, 7) =x(£)W(n), it may be demonstrated that with the product 
pairs FeymC€m, GEVmS€m, FeRmCem, GekmS€m, conditions (a), (b) cannot both be satisfied. 


Hence in (8.5)—(8.8) we put Fn =Gm=Fm=Gn=0. Thus the complete formal solu- 


tion is 


c(é, 7) = > * CmCem(E, 9)cem(n, 9) + CmCen(~, — g)cem(n, — 9) 


m=0 
- z: SmSe€m(E, 7)S@m(n, G) + SmS€m(E, — g)Sem(n, — 9). (8.9) 
m=1 


9. Boundary conditions. When the periphery of the plate is clamped, f =0f/d&=0 
at £=£). Thus for mode m corresponding to characteristic number a2, we have 
ConCeon(Eo, 7)Ce2n(n, 9) + ConCeon(to, — q)ceon(n, — g) = 9, (9.1) 
and 
ConCean(Eo, g)C€on(n, g) + F Kaa tis. — g)ceon(n, — q) = 0. (9.2) 
For Co,, Ce, to be non-zero, we must have 
(Ces, (£, q)Cean (Eo, gq) — Céon(fo, — g)Cean (Eo, q) |ceon(n, q)Céon(n, —g) = 90. (9.3) 


Hence the pulsatance equation is 


Ceéon(éo, —, )\Ceon (Eo, q) i Cen (Eo, wt g)Ceon(Eo, q) = 0, (n e 0), (9.4) 
or 
e .. — 
— [Ceon(E, q)/Ceon(é, — 9) lentes = 0. (9.5) 
dé 
This equation is satisfied when g=q;,, say, where s=1, 2, 3, - - - , and these roots 


define a system of confocal nodal ellipses. For the same values of g, by (9.3) we have 
the equations whose roots define two systems of nodal hyperbolae, namely, 


Céon(n, g) = 0, and Cés,(n, — gq) = 0 (n > 0). (9.6) 


For modes associated with the characteristic numbers den41, ben41, we take the respec- 


tive equations, with m=(2n-+1), 
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CmCem(€, q)cem(n, 9) + SmSem(~, — g)sem(n, — g) = 0, (9.7)* 
and 
SmS€m(E, )S@m(n, G) + CuCen(t, — q)cém(n, — q) = 0. (9.8)* 


Using the boundary conditions in (9.7) yields for a=4d2,,4;, the equations 


d 
[Cem(é, q)/Sem(&, = q) |e » ™ Q, (9.9) 
dé 
and 
cém(n, g) = 0, Sem(n, — q) = 0, (9.10) 


for g=4q;},., these being the roots of (9.9). 


Similarly for (9.8) we derive for a=be,41, the equations 


d 
= [Sem(E, g)/Cem(é, — 9) ]e-e, = 9, (9.11) 
ag 
and 
Sem(n, g) = 0, Cem(n, — g) = 0, (9.12) 


for g=4@,,,, these being the roots of (9.11). 

For characteristic number a@=)2,, the equations may be obtained from (9.5), 
(9.6) by writing Se for Ce, se for ce, G3n5 for G3.ns, m= 1. 

10. Elliptical ring plate. When an elliptical plate is clamped at its periphery and 
also at an internal confocal ellipse, it becomes a ring. Defining the inner and outer 
boundaries by &, and £9, the conditions to be satisfied are ¢ =0¢/dE&=0 at E=£, and £o, 
where 0S&<£ . If & were zero, the inner clamp would be on the interfocal line of 
length 2h. By Sec. 8, the formal solution for characteristic number dg,, is 


c@n) = [Cont Con (&, q) + Fo,Feyon(&, q) ]ceo,(n, q) 
+ [CoCeon(E, — ) + FonFekon(é, — q) |ceon(n, — 9). (10.1) 
Then at &=£ and at &=£,, we must have 
Ic Ceon(E, g) + FonFeyon(&, q) |ceon(n, q) 
+ [ Co,Ceon(E, cae q) + Fo, F ehon(é, = q) ]ceon(n, — q) = 0, (10.2) 


’ 7 Pa . . , 
[ConCean(, @) + FonFeven(, 9) |ceon(n, 9) 
ra] "i nh . , wa 
+ (C onC @on(~, — g) + FonFekon(&, — q) |ceon(n, —q)=0. (10.3) 
Thus there are four equations from which we can eliminate the four arbitrary con- 
stants, and so derive the pulsatance equation, and those for the nodal hyperbolae. 
Similar equations may be derived corresponding to @an41, Den41, Dene. 
11. Examples. 1°. Elliptical ring membrane. The first pulsatance equation is (6.3), 


* It should be observed that cézn4:(y, g) and Sézn4i(y, —g) have @=@e2n4:, while for séon4:(y, g) and 
Con4i1(n, —@), @=ben41. Similarly for Ce and Se. 
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and we have to find those values gm,» =Wy,,?p:/47 which satisfy it. In the absence of 
tabular values of the two functions involved, calculation of the values gq for the 
modes of lower order would be tedious. However, by aid of formulae asymptotic in 
k=+q'?, we can easily determine the approximate roots of (6.3). When g=k? is 
large enough? 

e | Kn cos 


eo v, — 6,), (11.1 
Fey q) ( ) 


(cosh pyle sin 


where K,, is a constant dependent upon g, v,=2k sinh &,, and 6,=(2m+1) tan-! 
(tanh 3£,). Using (11.1) in (6.3) leads to 


sin [(v) — v1) — (0 — 6:)] = 0, (73-2) 


(v9 — %1) = or + (09 — 6), (11.3) 
p integral 21. Thus (11.3) gives 
km, p = (pa + 00 — 01)/2(sinh & — sinh &:). (11.4) 


If e is the eccentricity of a confocal ellipse £ e-'=cosh & and, therefore, 
sinh &=e—1(1—e?)'/?, Substituting this into (11.4) yields 


9 


Qm,p = Rm,p 
222,32 2 22 f 2 2) 1/2 7 
= (pr + % — 61) e¢:/4 [eo + @1 — 2ee; — 2e1€9} (1 _ é:)(1 — €); l, (11.5) 


€o, €; being the eccentricities of the outer and inner bounding ellipses, respectively. 
The accuracy of (11.5) improves with increase in e and p. The second pulsatance 
equation (6.7) may be treated in a similar manner. 
2°. Elliptical plate. Referring to (9.5), we shall consider the modes of order 2n 
corresponding to the characteristic number d@z,. As before, owing to absence of tabular 
values, we shall assume that is large enough for (11.1) to be used. This formula 
applies when g>0, but we also require one for g<0. To derive this we write (+ 77) 
for £ in (11.1), since this substitution changes the sign of g in (3.3) of which (11.1) 
is an approximate solution when g is large and positive. As explained in the paper 
quoted in Footnote 3, it is also necessary to multiply by (—1)* and select the real part 
of the formula, since £ is real. Then if |£|>>47, we obtain 
Céon(E, — g) ~ Kenxi/(sinh £)"?, (11.6) 
where x,=cosh u cos ¢2,+sinh u sin den, u=2k cosh £, denm(4n+1) tan (tanh 36), 
and K2, is a constant dependent upon qg. Thus from (11.6) and (11.1), we get 
(=*) Céon(é, 9) (tanh £)'/? cos x (11.7) 
=) ———. S —__———_5 a 
Kon Ceon(E, oe q) X1 


with x =2k sinh £—(4n+1) tan-! (tanh 3£). Performing the differentiation indicated 
in (9.5), and equating to zero, leads after a little reduction to the pulsatance equation 


tan x = (1/x’ sinh 2&) — (x1’/xix’),  § = &o. (11.8) 
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If &) is large enough, x1’/x1x’~—1, while x’ sinh 2£>>1. Then (11.8) reduces to 
tanx~-1, or x=~(s — })z, (11.9) 


s being a positive integer such that x is large enough for the above approximate 
asymptotic formulae to be valid. Thus (11.9) entails 


2k sinh &) ~ (s — 4)a + (4n + 1) tan™! (tanh $65). (11.10) 
If £9 is such that 2k sinh £) ~ ke’, and tanh 3 ~ 1, we obtain 


dena = Bone & (s + n)?xe-2, (11.11) 
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A GENERALIZATION OF THE FINITE FOURIER TRANSFOR- 
MATION AND APPLICATIONS* 


BY 

IDA ROETTINGER 

Ann Arbor, Michigan 
Introduction. The purpose of this paper is to generalize and unify the methods 
used by Doetsch [8|t, Kniess [12], Koschmieder [13] and others to solve certain 
boundary value problems by finite Fourier integral transformations. To give an 
idea of the general method to be developed a formal solution of the following boundary 

value problem is here given by means of a particular transformation. 
Vite, 8) => YoAt, F+ x. O< 2 << e, O<— ¢ 
Y(x,0+) = 0, Yi(x,0+) = 0, 
Y(O +, t) = 0, Y a4 —,t) + AV (ar -, 8) = O, h 3, 

Let S} F(x) } =/GF(x) sin k,xdx =f,(kn), (compare section 1 below), where sin k,.x 
are the characteristic functions of y’’(x)+*y(x) =0, y(0) =0, y’(a) +hy(r) =0, that 
is k,, n=1, 2,---, are the roots of tan kr= —(k/h), k>0, then, if F(x), F’(x) are 


continuous and F’’(x) is sectionally continuous in (0, 7) 
S\F"(x){ = - ki fs| k,) + k,F (0) + sin kim |F’ (3) + hF (x) |, 
(see theorem 1 below). 
This last formula applied to the above problem with respect to x yields the follow- 


ing transformed problem: 


d*y.(Rn, t) 2 aa ) 
— = — kavi(kn, t) +S) 2;, 


dt? 
¥s(Rn, 0+) = 0, dy.(kn, 0 +)/dt = 0, where 
2 a 
S}x} = (sin k,w)k, — (xr cos kyw)k, . 


The solution of the transformed problem is 
9 
Val Rn, j= [(1 — cos k,f) ‘k, |S} xt. 
The inverse transformation S—'}f,(R,){ is given in terms of a Sturm-Liouville 


series (see section 2 below) as follows 
S14 f.(kn)} = >> N(k,)f(kn) sin knx in (0, 2), 
1 


where N(zk,) denotes the normalization factor of the characteristic functions sin k,.x, 
=, 2. :: 


* Received Oct. 18, 1946. Presented to the American Mathematical Society April 26, 1946. See also 
abstracts in Bull. Am. Math. Soc. vols. 50and 51. The author expresses her gratitude to Professor Ruel 
V. Churchill for his advice and encouragement in the preparation of this paper. The content of this paper 
is part of a dissertation submitted in partial fulfillment of the requirements for the degree of doctor of 
philosophy in the University of Michigan. 

+ The numbers in brackets refer to the bibliography. 
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Thus the formal solution of the above boundary value problem is given by 


20 
V(x, 2) = >> N(R,) [(sin kn) Rn — (x cos kaw)ka (1 — cos kt) sin k,x, 
1 

The type of problems to which the method is applicable is described in the sum- 
mary (see Sec. 11 below). 

Since there are at present no readily applicable general existence theorems for 
solutions of the class of boundary value problems considered here, the existence 
and uniqueness of the solution should be established in each particular case. However, 
this procedure is not carried through in the problems solved in this paper, since the 
main purpose here is to set up a method which leads quickly to a formal solution. The 
problems solved are chosen merely to give an illustration of the method. 

The method can be compared with that of the Laplace transformation (see re- 
mark in Sec. 3 below). As in the case of the Laplace transformation the present 
operational method does not claim to solve problems which cannot be solved by any 
other method. Its advantage lies, just as in the case of the Laplace transformation, 
in its direct, short and systematic approach. Problems in partial differential equa- 
tions which by a suitable change of variables are brought into a form in which vari- 
ables can be separated, can be solved directly. Tables of transforms, which are, except 
for normalization factors, tables of generalized Fourier coefficients, save time in the 
computation of solutions of practical problems. Furthermore this operational method 
has, due to certain theorems the advantage of systematically finding closed form solu- 
tions (see Secs. 6, 9,10, 11 below) and thus exhibiting qualitative properties of a solu- 
tion which may not readily be found by the usual methods. 

1. Definitions of the transformations S, C and 7. Definition 1. Let K={k} be 
a set of real numbers and Jet F(x) be a sectionally continuous function in (0, 7). The 
transformations S| F(x)} and C{ F(x)} are defined by the equations 


fe(k) (1) 


S\F(x) -f F(x) sin kxdx 
0 
C} F(x) = f F(x) cos kxdx = f.(k) (2) 


respectively. 

The transformations (1) and (2) are called the general finite Fourier sine and 
cosine transformation respectively, relative to the interval (0, 7) and the set K. The 
transformations (1) and (2) map the class of functions F(x) onto a certain class of 


functions f(R). 
A Sturm-Liouville problem. The characteristic functions of the following Sturm- 


Liouville problem are used in the definition of the transformation T below. 
y'"(x) + k®y(x) = 0, in (0,7), ) 
Li(y) = ary(0) + aey’(O) + asy(r) + aay'(r) = 0, (3) 
L2(y) b,y(0) + bey’(0) + dsy(r) + bay'(x) = 0, 


also with specializations on part of the constant coefficients a;, b;, i=1, 2, 3, 4. It is 
assumed that the L;, 7=1, 2, are linearly independent and that 
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aybe — dob = d3b4 — a4bz. (4) 


This last condition is to guarantee real characteristic values (see for instance [10] 
vol. 1 p. 352). 
The characteristic functions gx,(x), n=1, 2,---, of (3) are given by 


¢x, (x) = A sin kyx + Boos kyx, (5) 
where 


(A/B) 


ll 


— (a; + a3 cos knw — agk, sin kyr) /(Rad2 + a3 sin kaw + a4 ky, cos kyw) 
= — (b; + b3 cos kaw — dak, sin Rat) /(Rndb2 + 53 sin kaw + bk, cos kyr) 
and , are the roots of 
| (Rag + a3 sin ka + agk cos kor) (a; + a3 cos kw — agk sin kr) =0. © 


D(k) = | = 
| (kha + bg sin ke + byk cos kw)  (b; + bz cos kw — byk sin kr) 


By use of (4), Eq. (6) reduces to 
2k(arbe rl ob) = [ (asd; aa @b3) a k?(a4bo —_ a2b,) | sin kr 
+ [(asb1 — ayb4) + (a2b3 — a3b2) |k cos kr. (7) 


Remarks concerning D(k). (Compare [10] vol. 1, p. 362.) If D(k) =0 has no roots, 
then (3) has no solution. 

If and only if D(k) is of rank 0, that is, if every term in D(k) is 0, then the char- 
acteristic values are called double, since in this case (3) has two linearly independent 
characteristic functions with the same characteristic values, e.g. a, =b2= —a3—by = 1 
all other a;, 6; zero. In this case the characteristic values are 2, n=1, 2, - - - , and 
sin 2nx as well as cos 2nx are characteristic functions of (3). 

If D(k) is of rank greater than zero the characteristic values are single. This is the 
case in particular in the first, second and third boundary value problem. 

Definition 2. Let K= {kn}, n=1, 2, - - - , be the sequence of characteristic values 
of (3), that is, the roots of (7). Let 


gx, (x) = Ax, Sin knx + By, cos kyx (8) 
be given by (5), where 


A}; =-— ay + a3 COS kw - ask, sin k,@, 


— (Rade + az sin kaw + Agkp COS Ry). 


And let F(x) be a sectionally continuous function in (0, 7). The transformation 
T { F(x)} is defined by the equation 


T {F(x)} =f F(aden,(2)ax = f(kn). (9) 


The transformations S, C and T are linear. The interval (0, 7) is chosen for con- 
venience and without loss of generality. The transformation (9) maps the class of 
functions F(x) onto a class of sequences of numbers; except for a normalization factor, 
each sequence is the set of Sturm-Liouville coefficients of the corresponding F(x) in 
terms of the characteristic functions ¢g;,(x). The restriction to sectionally continuous 
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functions was made in view of applications. In general the functions F(x) need only 
belong to L*(0, 7): i.e. need only be Lebesgue square integrable. 

2. Inverse transformations. Since the function ¢;,(x) form a complete set of 
orthogonal functions on the interval (0, 7), the transformation (9) has an inverse 
transformation in the form of a Sturm-Liouville expansion. For, let F(x) be sec- 
tionally continuous and [F’(x) |? integrable in (0, 7) and at a point of discontinuity 
xo, let F(x) be defined as F(xo) =4[F(x0+0)+F(xo—0) |], O0<xo<z, and let N(k,) 
denote the normalization factor of the functions (8), then the expansion in the char- 
acteristic functions (8) of the function F(x) converges to the function F(x) in 
(QO, w); 1.e., 


3 


F(x) = Do N(Rn)f(Rn)¢e,(*), in (0, x), (10) 
1 
and the convergence is uniform and absolute in every closed subinterval of (0, 7) 
which does not contain a discontinuity (see [6] vol. 1, p. 371, compare also [5] pp. 
268, 272). Equation (10) gives a formula for the inverse T—'{f(k,)}, the function 
whose 7-transform is f(,). The inverse is unique. Thus 


T~{f(k,)} = SN (b)f(ka)on,(2), in (0,2), (11) 


where N(k,) is given by 


1/N(kn) = 3[e(a1 + a3) + aw(as + a4) — 0102 + aga,] 
+ sin kyw[k,w(a2a3 — aya4) + (1/ kp) (Redo, — 4,43) | 
+ (sin 2k, /4hq) [kn(as + a2) — (a; + a3)] 
+ w(k,d2d4 + a,a3) cos knw + (cos 2k,2/2)(aya2 — aay). 


Similarly the inverse transformations S{f,(k)} and C-{f.(k)} of the transfor- 
mations S{F(x)} and C{ F(x)} respectively are given in terms of a sine and cosine 
series respectively if K is the sequence of characteristic values of certain special 
cases of (3). e.g. a: +0, 630 all other a;, b;, zero, then K= {n}, n=1,2,---,and 
S{ F(x)} is the finite Fourier sine transformation as defined by Kniess (see [12]) and 
S—{f,(n)} is given by a Fourier sine series. 

3. Transformation of derivatives. The following two lemmas can be proved by 
integration by parts. 

Lemma 1. If F(x), F’(x), - - - F®*-»(x), are continuous and F®"(x) is sectionally 
continuous in (0, 7), r=1, 2,---, then 

a=r—l 
S{Fe@N(x)} = (— 1)"k?"S{ F(x) } + : > (— 1)*k*[k{Fer-2s 2)(0) 


s=0 


mae Fer -28 2) (qr) cos kr} a FCr-2s-1) (7) sin kr] (12) 


and 
s=r—1 


C{F@)(x)} = (— 1)rk*C{F(x)} + SS (— 1)8k**[RF@r-2-2)(4) sin ker 


s=0 


+ FOr-2s-) (77) cos kr —_ F@r-2e-1)(0) J, (13) 
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Lemma 2. If F(x), F’(x),--- F°-)(x), are continuous and F®-(x) is sec- 
tionally continuous in (0, 7), r=1, 2,---, then 
s=r—1 
‘IR Qr—1 ( ( r 32 af . " "3 sh2sh(2 2s—2 
S(Fe-)(x)} = (— 1)RC{ F(x) } + sin kw DS (— 1)*hF @r-2-2)(g) 
s=0 


1 


fe > (— 1)*R2 “1[Fer 28-1) (a) cos ka — F (27-2 (0) | (14) 


s=1 
and 
s=r—1 
C{F@r-)(x)} = (— 1)tk8 1S { F(x) } + sin ke DO (— 1)8t R221 PR Cr-2s-) (7) 
s=1 
s=r—l 
+ D> (= 1)%k®*[FC@r2-) (x) cos kw — F@r-2*-2)(0) J. (15) 
s=0 
Theorem 1. If F(x), F’(x), - - - F°*-?(x), are continuous and F®”(x) is sectionally 
continuous in (0,7),r=1, 2,--- ,andif w= %,(0), w2=¢%,(0), us=gux, (7), Us =X, (77), 


(g(x) as given by (8)), then there exist numbers \ and wp such that 
T F(x) } = (— 1)'k, T YF (x) } 
s=r—] 2 
+ >> (= 1)8h, [AL (FOr?) + pL a(R) |, (16) 
s=0 
where 
A = (tobe “++ U1b,) /(aybe — dob;) 
= k, + (a,b _— 9b) 1 (abe a a a4b1) k,, cos knw > (a3b, + k,,a4be) sin kn | 
and 
Bh = (303 + U4Q4) ‘(aybs a a2b}) 


- 
= (a,\b. — dob;) '[(ayay — G203)k, COS kaw + (a1a3 + k,a2d4) sin ku]. 


Remark. Property (16) exhibits the usefulness of the transformation T in obtain- 
ing quickly a formal solution of the following type of boundary value problem: 


m d?'F( x) 


>. Ae — =G(x), in (0,7), where 
r=0 dx?" 
G(x), the constants A», and the quantities L;(F°**-»), j7=1, 2, s=1, 2,---r-—1, 


are assigned. The constants a;, b;, 7=1, 2, 3, 4, are allowed to be different for each s. 
The transformation is also useful in boundary value problems in partial differential 
equations, where one or more variables behave like x in the above problem. Property 
(16) is analogous to the property of the Laplace transformation which expresses the 
transform of the nth derivative of a function in terms of the transform of the func- 
tion and the value of the Oth, - - - (n—1)th derivative of the function at 0 (see for 
instance [5] p. 8). 
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Proof of theorem 1. By combining (12) and (13) T{ F@?(x)} can be expressed as 


follows: 


s=r—l1 


T{F@(x)} = (- 1)"k, T{F(x)} + 3 (— 1)*h,, [uaF @r-2*-2) (0) 


s=0 
— UF Crm (Q) — ugh Or?) (qr) 4 ggF @r-22-)) (qx) J, 


If 
Fer-2s-2)(Q) = yw, FO@r-te-1)(9) = 


9 


Fer -28 2) (qr) = Ws, FOr-2s 1) (qr) = WwW, 


it remains to be shown that 


4 4 
UsW, — UjWe + UzW4 — UswW3 = X bw a;w; + n>, b; wi, 
1 


1 
which means it has to be shown that the following four equations are consistent: 
us = da, + why 
— Uy; = Ado + ube 
— Us = a3 + pubs 
Uz = Ad, + why. 


It can be left to the reader to obtain for instance \ from the first two equations and yu 
from the last two equations and to show that with these values of \ and yu the above 


equations are consistent. 
4. Illustration of the operational method. Particular cases of transformations. 
As an illustration of the use of the operator 7 a formal solution of the following 


boundary value problem is to be established. 
Vi2t+ Vy, = Q(x), O0<2x<z, 0< 4, 


V(x, y) is to remain bounded as y approaches infinity and V(x, 0+) =0, further- 
more, a¥(0+, vy) +0¥.(0+, y) =£, where E is a constant, and c¥(r—, y) +dV,.(r—-, y) 
=0, where a, b, c, d are given non-zero constants such that ad —bc 0. 

This problem calls for the operator (9), where ¢gy,,(x) are the characteristic func- 
tions of y’’(x) +ky(x) =0, ay(0) +by’(0) =0, cy(r) +dy'(r) =0. From (7) follows that 
the characteristic values k, are the roots of tan kr=k(bc—ad)/{ac+k*bd), k>0O. 


From (9) it follows that 
T{F(x)} = f "F(x)(a sin k,x — bk, cos knx)dx = f(Rn), 
n=1,2,---,and (11) yields 
T-"{ f(kn)} = > N(kn)f(Rn)(a sin k,x — bk, cos kax), in (0, x), (17) 
1 


where N(k,) is given by 


1/N(R,) = 3(ra? + rhb? — ab) + [(hab? — a’)/4k,| sin 2kaw + (ab/2) cos 2kpm. 
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Equation (16) with r=1 yields 


T{F’'(x)} = — kaf(kn) + kn[aF(0) + bF’(0)] 
— [kab/(c sin kaw + Rnd cos kam) ][cF(x) + dF'(x)]. (18) 


Formula (18) applied to the above problem with respect to x leads to the follow- 
ing transformed problem: 


Vuln ¥) — Rab(Rny ¥) + BnE — q(Rn) = 0, 


W(Rn, y) remains bounded as yo and y(k,, 0+) =0, where ¥ and g stand for the 
T-transforms of V and Q respectively. 

The solution of the transformed problem, that is the coefficient in the Sturm- 
Liouville expansion of the solution is 


(kn, y) = [1 — exp (— kay) ] [An — (hn) ]/he 


and according to (17) the formal solution of the given problem is 
U(x, y) = D2 N(kn){ [1 — exp (— ny) ][RnE — q(n) ]/kn} (a sin kn — bhp COS ky). 
1 


In the following sections a special case of the transformations S and C will be 
needed (for additional special cases see section 12). The kernels sin kx and cos kx 
respectively of the transformations are the characteristic functions of y’’ (x) +k*y(x) =0 
y(0) =0, y’(r) =0; and y’(0) =0, y(7) =0 respectively. 


In these two cases K = {n— , ,n=1,2,---. The transformations are 
{ F(x) } -{ F(x) sin (n — $)xdx = f,(m — 4), n=1,2,--- (19) 
0 
C{F(x)} -{ F(x) cos (n — 3)xdx = f.(n — 3), 2=i1,2,---. (20) 
0 


The inverse transformations are given by 


as { fa(n _ 1) } = 


a | 


(n — 3) cos(n—4)x in (0,2). (22) 


Thilo — 3)sin(n—4)x in (0, 7) (21) 
=D id 
1 


C-1{ f.(n — 4)} 


And for the transform of an even derivative (12) and (13) yield 


S{FOC”(x)} = — (n — 4)**f,(n — 4) 
s=r—l 
$°D (= In — P"[(— YHFO* Oe) + (n= PFO] (23) 
s=0 
C{F@(x)} = — (n — 4)*f.(n — 3) 
s=r—1 


+ DO (= 1) — )™[FO-*-V) + (nw = A)(— 1) FO) J. (24) 


s=0 
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Remark. In the above listed cases D(k) is of rank greater than zero. If the rank of 
D(k) is zero it merely has to be kept in mind that the sine as well as the cosine coeffi- 
cients in the Sturm-Liouville expansion of the solution have to be found. An illus- 
ration of the operational method in this case is given in section 10. 

5. Some properties of the transformations S and C when K={n-—} 
n=1, 2,---+.Theorems 2 to 5 below are analogous to theorems which holds for the 
finite Fourier sine and cosine transformations (compare [5] sections 95 and 96 also 
[19] chap. II sect. 2.1). and can be proved in a similar way. 

Definition 3. Let the function F(x) be defined in (0, »). By the odd antiperiodic 
extension F,(x) of F(x) with antiperiod 2p is meant Fi(x) = F(x), in (0, p); Fi(— x) 
= — F(x), Fi(x+2p) = — F(x) in (— @, &), 

Definition 4. Let the function F(x) be defined in (0, ). By the even antiperiodic 
extension F,(x) of F(x) with antiperiod 2p is meant F(x) = F(x), in (0, p); F2(—x) 
= F,(x), F2(x+2p) = — F(x), in (— ©, &). 

Theorem 2. If F(x) is sectionally continuous in (0, 7) and if Fi(x) is the odd anti- 
periodic extension of F(x) with antiperiod 27 and if @ is any constant, then 


sin (n — 4)aS{F(x)} = 4C{Fi(x + a) — Fi(x — a) , 
Theorem 3. Under the same assumption as in theorem 2 
cos (n — })aS{F(x)} == AS{Fi(x + a)+Fi(x — a)}. 
Theorem 4. If F(x) is sectionally continuous in (0, 7) and if F2(x) is the even anti- 
periodic extension of F(x) with antiperiod 27 and if a is any constant, then 
sin (n — 3)aC{F(x)} = 3S{Fo(x — a) — Fo(x + a)}. 
Theorem 5. Under the same assumption as in theorem 4 
cos (n — 3)aC{F(x)} = 4C{F.(x — a) + Fo(x + a)}. 


Similar theorems in the case when K is the sequence of characteristic values of 
(3) and for more general K can be proved by the use of almost periodic functions (see 
[15]). 

6. Vibration of a horizontal string with one end fixed and one end sliding. 
Compare [5] sections 98, 99.) Let the end x=0 of a string be fixed and the end x= 
be looped about a vertical support along the line x =z. If a constant upward force acts 
on the loop, the displacements Y(x, t) as the string is released from the position Y=0 
satisfy the conditions: 


Vi(x, t) = a®V,,(x, t) + g, 0O<x<z, 0<4, 
Y(x,0+)=0 Y,(x,0+) =0, 
Y(0 +, 4) 0 Y (x —, 2) — b, 


where g is the acceleration of gravity and b is the magnitude of the vertical force 


divided by the tension. 
The S-transformation (19) applied with respect to x gives, according to (23), the 


following transformed problem: 
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dy(n — 3, 2) 
— — + (n — 4)*a*y(n — 3, 2) + (— 1)"*10% — gS}1; =0 


dt? 
y(n — 3,0+) =0, y(n — 4,0+) =0. 


The solution of this transformed problem is 
1 


y(n — 3, t) = [1 — cos (nm — 4)at]} [g/a2(n — 4)3] + [b(— 1)"/(n — 4)?]}. 


From the table in section 12 below it follows that 


S-'{ [b(— 1)"/(n — 4)2] + [¢/a2(n — 3)*]} = — bx + [gx(2e — x)/2a?] = F(x). 
Using theorem 3 the solution of the given problem can be written in the form, 
Y(x, t) = F(x) —43[Fi(x+at) + Fi(x —at) J, where F,; denotes the odd antiperiodic 


extension with antiperiod 27 of F(x). 

Similarly the problem with the end x=0 sliding and the end x=7 fixed can be 
solved by use of the C-transformation (20). 

7. The convolution in the case K = }n—}}. The purpose of this section is to give 
formulas for the product of two transforms in terms of one transform. Kniess (see 
[12], compare also [19] 2.1) and Doetsch (see [8]) give such results for the finite 
Fourier sine and cosine transformations, i.e. when K = {n , n=1,2,---+-+. Here the 
analogous results are given for the S- and C-transformations when K=j{n—}} 
(compare also Koschmieder [13]). Analogous results when K is the set of char- 
acteristic values of certain Sturm-Liouville problems can be obtained by the use of 
almost period functions (see [15]). The proofs of theorems 6-9 are analogous to those 
of Kniess. 

Definition 5. If F(x) in (—2z2, 27) and G(x) in (—7, 7m), are sectionally continu- 


ous, then the function 


F(x) *G(x) -{ F(x — y)G(y)dy 


is called the convolution of F and G on the interval (—7z, 7) (compare [12] p. 270 and 


[5] p. 274.) 
Lemma 3. If F(x) and G(x) are sectionally continuous and F(x+27) = — F(x) and 


G(x+27) = —G(x), then 
F(x) *G(x) = G(x) +F(x). 
Theorem 6. If F(x) and G(x) are sectionally continuous and even functions and if 
F(x+27r) = — F(x), then 
C{F(x)}C{G(x)} = 4C{ F(x) +G(x)}. 
Theorem 7. If F(x) and G(x) are sectionally continuous and ‘odd functions and if 
F(x+27) = — F(x), then 
S{F(x)}S{G(x)} = — $C) F(x) «G(x) I 
Theorem 8. If F(x) is an even and G(x) is an odd sectionally continuous func- 
tion and if F(x+27) = — F(x), then 
C{F(x)}S{G(x)} = 3S {F(x) +G(x) }. 
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Theorem 9. If F(x) is an odd and G(x) is an even sectionally continuous function 
and if F(x+27) = — F(x), then 


S{F(x)}C{G(x)} = 4S{F(x) +G(x)}. 


Remark. If in the above four theorems G(x) satisfies the condition G(x+27) 
= —G(x), then according to lemma 3 the convolution is commutative. That being the 
case theorems 8 and 9 say the same thing. 

Example. Given q.(n—}) =(n—}3)[1—(n—4)?]}-, find Q(x). ge(n—4) =S{x}X 
S}sin x; (see tables section 12). According to theorem 7 if F(x) =x in (0, 7) and 


F(x + 2x) = — F(x) \ 
‘ In (— oO, 2 ) 
G(x) = sin x 
1 z-f : rT ; 
Q(x) = - : lf (24 — x + y) sin ydy +f (x — y) sin vay 
=mr—x+sinx. 


8. Basic problem. In section 9 an application of the solution of the following 
problem will be made. 


0} 


4 m geny 
+ >} (-1)"——- = 0, 0<2x<z, 
at oni 


— o<st 
Ox" 


g2m—2e-2)V (0 +, 2) gem—2e-DV(e —, 2) 


: $$ _______—— 0, 
Ox2m—28-2 2m—! 
s=0,1,2,---,m—1, Y(x,0+) = 1. 


The S-transformation (19) with respect to x gives according to (23) the following 
transformed problem: 


m 


dy(n — 3, 2) 
— "7 "+ y(n — 3,2) > (n — 3)* =0 
at 


r=] 
y(n — 3,0+) =S{1} = 1/(n — }). 


The solution of this transformed problem is 


y(n — 3, t) = (mn — 2) exp | - DY un - p*|. 
r=1 


According to (21) the solution of the basic problem is thus 


V(x, 2) = (2/r) >> (n — 4)— exp | - a t(n — 3 *| sin (m — $)x. 
n=1 r=] 


9. A generalized heat conduction problem. A formal solution of the following 


generalized heat conduction problem is obtained by use of the operator C when 
K={n—3} (see 20)). 
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_ aU m g2nU . 
L(U) =— + Ao(t)U(x, 2) + D> Aa,(t)(— 1)" = Q(x, 2), 
ot gust ox?" 


0O<x<z, i> 0, A2,(t) > 0, 


g@m—22-1) (0 +, 2) 
_- ———-= o2 Bom 2s—1(), 











Ox2 m—2s-1 
s=0,1,2,---,m-— 1, 
Q(2m—28-2) [7 (ar -, t) 

————- = Daw o5_9(t), 


ra) x2 m—2s—2 





U(x«x,0+) = F(z). 


This problem can be resolved into 2m+2 problems each of which has 2m+1 homo- 
geneous conditions and one non-homogeneous condition. The sum of the solutions 
of these problems is the solution of the given problem. There are four essentially dif- 
ferent types of problems, a formal solution of which is given in the following. 


Problem I. L(U) =Q(x, t), 0<x<7, 0<t, 0<A2,(t), 
g(2m 2e- NU (0 +, t) gm 28 U(x _v t) . 
SO = ), ad 


4 -2m—2s—1 ) -2m—28—2 : 
Ox Ox 


s=0,1,2,---,m—1, U(x,0+) = 0. 


The C-transformation (20) applied with respect to x leads according to (30) to 
the following transformed problem: 


du,(n — 3, t) - 
—_——— + | 40 + >) A2(t)(n — »*| u(n — 3, t) — gen — 3, t) = 0, 


r=1 
u(n—4,0+) =0. 
This equation is of the type y’(x)+P(x)y(x)—Q(x) =0 with initial condition 


y(0) =0. The solution thereof can be written in the form 


y(x) -f exp [— p(x, w)|Q(w)dw, where p(x, w) -{ P(v)dv. 
0 w 


Hence the solution of the transformed problem can be written in the form 


m 


u(n — 3,t) = f ‘* [— Ao(t, v) | exp | - >> A2,-(t, »)(n — *] ge(n — 3, v)d, 
0 r=1 
where 
Ao,(t, v) =f AsGvas. 
Now 
exp | - » A>,(t, v)(n — a* | = (n — A)S{Y (x, A2,(t, v))}, 
r=1 


where Y is the solution of the basic problem of section 8. And according to (15) with 


k=n—3 and r=1. 





1947] FINITE FOURIER TRANSFORMATION 309 


(n — 3)S\Y(x, Aor(t, v))} = CLY (x, Aae(t, v))}, 


so that 
u(n—4,i) = foe [— Ao(é, ») Ic{ Y .(x, A2,(t, v)) 1C{Q(x, v) \ dv. 
Using theorem 6 the solution of problem I can be written in the following form: 
U(x, t) = f ‘as [— Ao(#, v) ][4¥ -(x, As-(t, v)) *Qo(x, v) ]dv, 
0 


where 
Qo(x, v) = Q(x, v) 0 
Qo(— x, ») = Qo(x, 2), —fT 

Problem II. L(U)=0, 0<x<7, 0<#, A2,-(t) >0, 
gc2m-2- (0 +, 2) 


pape’ Bom—os—-1(1) for s = i, 
C x m—<8— 


= 0, for s ¥ 1, 


gm -28 DU (x —, t) 





0, U(x,0+) =0, 


s=0,1,2,--+,m—1. 


According to (24) the C-transformation applied with respect to x leads to the fol- 
lowing transformed problem. 


m 


du(n — 3, t) : 
ee | 40 + D A2-(t)(n — 3 «| ue(n — 3, t) 
( r=1 


+ ¥ As (i)(— *(— D4(n — )*Barsea(t) = 0, 


r=0 
u(n — 4,0+) = 0. 


With the notation of problem I the solution of this transformed problem reads: 


u(n — 3,t) = f ‘a [— Ao(t, ») ] exp| - d Ant v)(n — »*| (n — 3)*(— 1)* 
0 r=1 


- YS (= 1)*Bar—2i-1(0) A 2(0) do. 
r=0 
By use of the solution Y(x, ¢) of the basic problem and (15) with r=i+1 and k=n—} 


(n — 3) exp | - > Axl, o)(n — *| = (n — )*S{Y(x, An(t, »))} 





acetDV (x, A2,(t, ” 


Ox2tt1 


=(- nec 
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Hence the solution of problem II reads 


t i HOY (x, A x(t, v)) m 
U(x, t) -{ exp [— Ao(t, eT (x (— 1)"Ao,(v) Bo,_2; ()) de 


grit r=0 


Problem 1II. L(U) =0, 0<x <7, 0<t, A2,(t) >0, 


Q (2m—2s-1) 17 (0 +, #) 
- - 0, s=@, L.2,°-: ,m=—1, 


0 x m—2s—1 


g(em 2s DU (9 —, #) . 
; ee ake = Dom os2(t), S$ = 1, 
Ox- m—28-—-< 


0 $s ¥ 14, 


(x, 0 +) 0, 


According to (24) the transformed problem reads: 


du(n — 3, t) 
cei aie + | Aut + a t2,(t)(n — 3) ax — 4,8 
at r=1 


= 0 


+ S) (= 1)"A2,(é)(— 1) — 1)"(m — 4)? Do, oi2() 
Tr=VU 
u(n — 3,0+) = 0. 


The solution thereof in the notation of problem I is 


9541 


t 
u.(n — 3, t) -f exp [— Ao(t, »)] exp | - 2a Ao,(t, v)(n — 3)? | 1)"*+*(n — 3)? 
. b Ds (- 1)*A o,-(v) Do, 2; o(v)dv. 


r=(0 


By use of the solution Y(x, ¢) of the basic problem and (23) with r=7+1 


(n — 4 exp | - 2: Ao,(t, v)(n — 3)? ‘|= (n — 3)+°S{ V(x, Aar(t, 2)) } 


collie (x, Aor(t, v))) 
= (— 1) +19 — ~ 3 f° 


ox” 


Using theorem 3 
gee of ee Avis 2))\ 


i 
eee sf 


Ox? 


so that the solution of problem III reads 
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U(x, t) [ ‘aa [— Ao(t, v) | ili ssc ” | > (- 1)"Anl0) Davis) | dv, 
J 9 dy" r=0 
yor. 
Problem IV. L(U) =0, 0<x<z, 0<t, Ao,(t) >0, 
g@n=-2e-DU (0 +, 8) Geo U(e'—, 6 
= a ar ia 
s=0,1,2,---,m—1, U(x, 0+) = F(x). 
According to (24) the transformed problem reads 
1 m 
an a8) + EG + 2, Aar(t)(m - »*| u(n — 3,t) = 0, 
u(n — 3,0+) = f.(n — 3). 


The solution thereof in the notation of problem I can be written as 


| 
II 


fn — 4) exp [— Ao(t, 0)] exp | - > A2,(t, 0)(n — *| 
r=1 





-_ AV (x, Ao,(t, 0)) 
= exp [— Ao(t, 0) |C} F(x) cf ~ ; 
O° 


since 


(n — 1)S} Y(x, Ao,(t, 0)) 


C jue. A o(t, a 
Ox 
according to (15) with r=1 and k=n—}. 


By theorem 6 the solution of problem IV can be written as 


1 AY (x, A2,(t, 0)) ] 


exp | - > A2,(t, 0)(n — »*| 


r=1 





eee Wel 
2 Ox o(z) 


U(x, t) = exp [— Ao(#, 0) | 


where F(x) = F(x) in (0, 7) and Fo(—x) = Fo(x) in (—7, 7). 

It can be seen now that the solution of each problem is expressed in closed form in 
terms of the solution of the basic problem of section 8. Thus the solution of the given 
problem being the sum of the 2m+2 problems is expressed in closed form in terms of 
the solution of the basic problem. 

10. A problem in heat conduction. As a particular case (m=1) of section 9 the 
solution of the following problem can be obtained. 

To find the temperature distribution U(x, ¢) in a slab of length 7 with a heat 
source inside, The end x=7 is kept at a temperature D(t) and radiation through the 
end x =0 at a rate B(t} takes place. Furthermore there is radiation through the lateral 
surface and the thermal diffusivity K(¢) depends on time. The initial temperature of 
the slab is F(x) .(Compare also [1].) 
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The mathematical formulation of this problem is the problem of section 9 with 
m=1, A2(t)=K(t), Ao(t)=A(t), Bi(t) = Bt), Do(t) =D(8), ie. 


oU eU 
— -— K(é) A(t)U(x, t) = Q(«, 4), 0 <= <s¢, 0 < ¢ 
Ot Ox 
dU(0 +, 24) 
Ox 





= Bi) U(r-,t) = Di, 
U(x,0 +) = F(x). 
Put 

_— t “— t 

K(i, v) -{ K(x)dx, A(t, v) -f A(x)dzx. 
The solution of the basic problem of section 8 with m=1 is 

V(x, #) = (2/m) D (n — §)~ exp [— (n — )%] sin (n — Ba. 
1 


The solution of the heat conduction problem can be written in closed form using 


Jacobi’s 2-function. 


For 
oY (x, 2) 





= (2/m) >> exp |— (n — 4)*t] cos (n — 4)x 
x 1 


= (1/r)82(x/2, exp [— t]), (see [12] p. 464.) 


Thus the solution reads 


U(x, t) = f P= [— A(t, v) ][(1/22)82(x/2, exp {— K(t, v)})*Qo(x, v) |dv 
~ (1m) f po [— A(t, v) }82(x/2, exp { — K(t, v)})K(v)B(v)dv 


t oe Ove a 
~ (1/x) f exp [— A(t, 0)] — (y/2, exp {— Ke, 9)})K@)D)dn 


y=r-—% 
+ exp [— A(t, 0) ][(1/2x)82(x/2, exp { — K(t, 0)}) «Fo(x)]. 
Remarks. Since the 3-functions are tabulated the above form of the solution can 
be used to determine by mechanical integration the numerical values of U(x, #) for 


given values x and ¢. 
For questions on uniqueness and existence of solutions in heat conduction prob- 


lems the reader is referred to [9]. 
11. An illustration of the method in the case of double characteristic values. 


Summary of the operational method. A formal solution of the following boundary 
value problem is to be found. 
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L(U) = U(x, t) — K()U 22(x, t) + ADU (x, t) = Q(x, d), 
0O<x<7, 0<4%, A(t) > 0, K(t) > 0, 
L,U) = U{0+,) + Ur —, fd) = Bd, 
LU) = U(0+, t) + U(x -, t) = DY), 
U(x,0 +) = F(x). 

As kernel of the transformation (9) to be applied to this problem the character- 
istic functions of y’’(x)+k?y(x) =0, y’(0)+y’(r) =0, y(0)+y(7r) =0, are used. In 
this case the determinant D(k) is of rank zero. The double characteristic values are 
2n—1, n=1, 2,---,and sin (2n—1)x as well as cos (2n—1)x are characteristic 
functions. This means the sine as well as the cosine coefficients of the solution of the 
above problem have to be found. The formulae to be used are: 


S{F’'(x)} = — (2n — 1)*f,(2n — 1) + (2n — 1) [F(O) + F(w)], (25) 

C{F’(x)} = — (2n — 1)*f-(2n — 1) — [F'(0) + F(m)], (26) 

C{F'(x)} = (Qn — 1)f,(2n — 1) — [F(r) + F(O)]. (27) 
The last three equations follow from (12), (13) and (15) respectively with r=1 and 
k=2n—1, n=1, 2,--.-. Relation (11) yields 


S“{f(2n — 1)} = (2/n) O fal2n —1)sin(2n—1)x in (0,2). (28) 


Furthermore the solution of the following auxiliary problem will be useful. 
Y.(x, t) — Y.:(x,i) = 0, in (0,7), Y(0+,4+ Yr -,#) =0, 


V(O+, t)+V(r-, t)=0, V(x, 0+)=1. From (25) it follows that [dy,(2n—1, t)/dt] 
+(2n—1)*y,(2n—1, t)=0, y.(2n—1, 0+)=S{1 } =2/(2n—1). Hence y,(2n—1, #) 
=2 exp |—/(2n—1)?](2n—1)—, which gives the sine coefficient of the solution of the 
auxiliary problem. The cosine coefficient is zero since C{1} =0 when k=2n—1; so 
that the solution of the auxiliary problem is according to (28) 


Y(x, #) = Q/n> 2(2n — 1)! exp [— (2n — 1)*t] sin (2n — 1)z, (29) 


VY (x, t) = (2/m)82(x, exp [— 4#]), see [18] p. 464. (30) 


The following formulae by Kniess (see [12]) will also be used. If K= {mn} then a) 
if F(x) is even and periodic with period 27 and G(x) is odd, then 


C{F(x)}S{G(x)} = 4S{F(x) +G()}; (31) 
b) if F(x) and G(x) are even and F(x) is periodic with period 27, then 
C{F(2) }C{G(x)} = AC {F(x) +G(a)}. (32) 


The solution of the given problem can be written as the sum of the solutions of 
the following four problems, where 0<x<7, 0<#, A(t)>0, K(#)>0. 
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I. L(U) = Q(x, 2), L,(U) = LU), U(x,0+) = 0. 
Il. L(V) = LU), U(«*,0+) = 0, L,(U) = B(t). 
III. L(V) =L,(U), U(x,0+)=0, L(U) = Did). 
IV. L(V) = L,(U) = L(V) = 0, U(x,0+) = F(x). 
The sine and the cosine transforms using (25) and (26) respectively of the solution 
of each of these problems have to be found. The sine transform of the solution of II 


and the cosine transform of the solution of III are zero. 
With the notation of section 10 it follows from (25) that the sine transform of the 


solution of I is 


t 
u,(2n — 1, t) -{ exp [— A(t, v) | exp [— (2m — 1)?K(t, v) ]q.(2n — 1, v)dv. 
0 
Relations (29) and (27) yield 
exp [— (2m — 1)*t] = [(2n — 1)/2|S{Y(x, )} = 4C{Y.(x, dD}. 


Using (31) and (30) the sine part of the solution of I can be written as 
t 
(1 jan) f exp [— A(t, v) ][#2(x, exp (— 4K (tz, v)})*Qo(x 2) |dz, 
0 


where Q is extended so as to be odd with respect to x. 
From (26), (29) and (27) it follows that 


1 ; = ar 
u-(2n — 3, t) = —{ exp [— A(t, v) ]JC{V.(x K(t, v)){C{Q(x 2») fade. 


( 


According to (32) and (30) the cosine part of the solution of I can be written as 
t 
(1 2m) f exp [— A(t, v) ][#e(x, exp | — 4K(é, ») }) «Q.(, 2) ]do, 
0 


where Q is extended so as to be even with respect to x. 

It may be left to the reader to show that in a similar way the sine and cosine 
parts of the problems II-IV can be expressed in closed form in terms of J2 by using 
relations (25) to (32). Thus the solution of the given problem can be expressed in 
closed form in terms of Jacobi’s #2 function as follows: 


t 
U(x, t) = (1/2x) f exp [— A(t, v) ][do(x, exp | — 4K(i, v) }) *Qo(x, v) ]do 


0 
t 
+ (1/2n) f exp [— A(t, v) |[do(x, exp | — 4K(t, v)}) «O(a, v) ]do 
0 
, —_— _—_ 
_ (1/n) f exp [— A(t, v) ]82(x, exp | — 4K(t, v) }) B(v) K(v) do 
0 
P i _ 
= (1/m) f exp [— A(t, v) ]82(x, exp | — 4K(t, v)})D(v) K(v)dv 
0 


+ (1/21) exp [— A(t, 0)][d.(x, exp | — 4K(t, 0) })*Fo(x) | 
+ (1/2) exp [—A(t, 0) ][8o(x, exp | — 4K(t, 0)})*F.(x)]. 
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Summary of the operational method. The method is applicable to the following 
type of boundary value problem: 


y (— 1A» ee = Ge) in (0,5), 


r=0 
where G(x), the quantities 
a 2r—28 2)(0) + aoF @r “2s—1)(() + a3F (2r-28- 2) (ar) o a 4F (2r-—2s-— 2) (a) 


and 
b.F 2r-28 2)(0) + boF er 28 (0) + b3F (2-2-2) (ar) +. b4F (27—2#-2) (qr) 


and the constants A», are assigned, s=0, 1, 2,---,r7r-—1. The constants a,;, };, 

i=1, 2, 3, 4 are allowed to be different for each s. The transformation is also useful 

in boundary value problems in partial differential equations, where one or more 

variables behave like x in the above problem. In order to obtain a formal solution 

quickly, and possibly a closed form solution, the following procedure is recommended : 
1° Set up (3) corresponding to the given problem. 

2° Find the roots of D(k) =0 (see (7)). 

3° Find the rank of the determinant D(k). 

4° If the rank is greater than zero set up T{ F} using (9). If the rank is zero set 
up 7:| F} and T:{ F}, where 7; aad 7; have the two independent characteristic func- 
tions of (3) as respective kernels. 

5° Find T} F?"} using (16). 

6° Apply T} F?"} to the boundary value problem (which may be resolved into 
several problems each of which has only one non-homogeneous condition) and find 
f(R,), the transform of the solution. 

7° Try to obtain the inverse T-'{f(k,)} in closed form by application of the theo- 
rems of sections 3, 5, 7 or by use of tables of transforms or a combination of both. 

8° If 7° does not lead to the solution, use (11) to find T-"{f(k,)} and obtain the 
solution in series form. 

Remark. Formulae for a few special cases of T can be found in the tables section 
12. 

12. Tables. Tables A—D contain a few examples of transforms. It would be de- 
sirable to have extensive tables of transforms, since they would help in obtaining 
closed form solutions of boundary value problems (see the example of section 6). 

Table E below contains a list of transformations which are special cases of the 
transformation T (see (9)) of section 1. In each case the formula for the transform of 
an even derivative (see (16)) as well as the inverse transformation (see (11)) are given. 
The completion of this list of special transformations is left to the reader. 





~ 
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, A. Tables of S-transforms. 
ee (ee Salk) =S{ F(x)} F(x) 
:. (1—cos kr)k7 1 
i k sin kr—k™"n cos kr x 
gf 
3. | [(—1)"(2r)(1—cos kwr)/k*1]+cos ke D> (—1)rt*[(2r) tx#2/ (254-2) lk | a 
| sal 
r—l 
| +sin ka >> (—1)rtett[(2r) fr2*t1/(2s-+1) le] 
s=0 
| r—l | 
4. | sin ke >, (—1)rte+1[92(27—1) !/(2s) 1k2*] ar 
+cos kw >> (—1)*t*[n?"+1(2r—1) 1/(2s-+1) 1h | 
s=0 
3. (k?—a?)—“(a cos aw sin kr—k sin az cos kr), a#k sin ax 
6. (k?—a?)“(k—k cos kw cos aw—a sin kr sin ar), a#k cos ax 
B. Tables of C-transforms. 
Pee — fhk=CiP@)} = FR) > 
1. | k*sin kr 1 
2. | ko wsin kx—(1—cos kr)k? x 
3. | sin kw >, (—1)rt*[w2*(2r) !/(2s) tk? | xr 
rl 
+cos kr), (—1)rt#+1 [x +1(2r) !/(2s+1) lk-**] 
sal) | 
r—l | 
4. | [(—1)"(2r—1) !/k*"]+-cos kr D> (—1)rt#+1[(2r—1) !x/(2s—2) 1k] | x21 
s=0 
ame | 
+sin kr >, (—1)rt#+1[(2r—1) !%*1/(2s—1) 1k? 2 | 
s=0 
2. (k?—a?)“1(—a+k sin ka sin ar+a cos kx cos ar), ak sin ax 
6. | (#—a?)“(k cos am sin kr—a sin ax cos kr), a¥k cos ax 
C. Table of S-transforms when K={n—}},n=1,2,-+-. 
f(n—4) =S{ F(x)} F(x) 
[. (n—})7 1 
2. | (—1)"(n—4) x 
a (n—})-3 x(2xr—x)/2 
1—cos x 


4, (n—3)71(1—(n—})?) 
5. | (n—4)[(n—})?—a? 7 


6. (—1)"*1a[a*+(n—4)?]71 


cos a(r—x)/cos ar 


sh ax/cosh ar 


D. Tables of C-transforms when K =\n—}},n=1,2,+-> 


f(n—4) =C{ F(a) } F(x) 
| (—1)"4#(n—4) Ser | -+s)0°”~CU ae ae 
(—1)"*4x(n—})1—(n—})? x 
aX 


| . 
(n—4)7[1-—(n—4)*]71 w—x—sin x 


[(n—4)*—a?]7 
6. (—1)"41(n—4) [(n—})?—a?]71 
ds [a?-+(n—4)?)> 
8. (—1)"**(m—4) [a?+(n—})?}7 


1 

2 

S:, | (n—})? 

4 

5 sin a(r—x)/cos ar 
cos ax/cos ar 

| sh a(r—x)/cosh ar 


| cosh ax/cosh ax 
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E. Table of some special transformations. 


S{F(x)} -f F(x) sin nxdx = f,(n), n=1,2,--- (33) 
0 


S-1{ f,(n)} = (2/n) > flr) sinnx in (0, 7) 


s=r—l 
S{FOM(a)} = (= nfm) + YO (= 1m [Rer-e-9(O) — (— 1) Fer-2-9(q) 
s=0 
C{F(x)} = f F(x) cos nxdx = f.(n), 2 =0,1,2,---> (34) 
0 


C-{ f.(n)} = (1/)f-(0) + (2/r) > f.(n) cosnx in (0, x) 
1 


s=r—1 


C{PO(x)} = (= nfm) + YO (= 1)n™[(= 1) Fert) (g) — Fere-0)] 
s=0 


S\F(x)} = f F(x) sin (n — $)xdx = f,(n — 3), n=1,2,--- (35) 


S-1} f.(n — 4)} (2/r) >> f.(n — 4) sin (n — 4)x in (0, x) 
1 
S{FO”(x)} = (— (n — 3)*)"f.(n — 3) 


+ DO (= 1) — He[(— Fe-0@) + (n — HREM“) 


s=0 


ChF(x)} -f F(x) cos (n — 4)xdx = fn —}), n=1,2,--- (36) 
0 


C-'} f.(n — 4)} = (2/n) > f.(n — 4) cos (n — 4)x in (0, x) 
1 


CiFC(x)§ = (— (n — 3)%)"f-(m — 3) 


s=r-l1 


+4 >» (- 1)*4(n a 4) 2e [Fer 2s- 1) (()) + (n — 1)(— 1) "F (2r- 28-2) (qr) | 


S{F(x)} = f F(x) sin k,xdx = f.(Rn), gmi1,2,-+-, (37) 
0 


where &, are the roots of tank = —(k/h), k>0, h¥0 


oo 


> N(kn)fe(Rn) Sin knx in (0,7) and 


1 


(3/2) — (4k,)— sin 27k, 


Sf f,(Rn) } 


1/N(k,) 


Il 
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S{Fe(x)} = (- — 


s=r—l 


+> (- [sin kw [FO2r-2*-) (ar) + AF Or-2*-2) (7) | + hy r-2*-2)(0) | 
s=0 
S{F(x)} wile F(x) sin kn(m — x)dx = f.(Rn), n=1,2,--- (38) 
0 
5-1} .( m2 N(kn)fe(Rn) sin ka(w — x), in (0,7), 
1 


k, and N(k,) as in (37) 


S{ FO? (x)} = (— ka) fel Rn) 


s=r—1 
+d (—1)*H2, | sin kyw[F(7-2*-) (0) — hF @r-2#—-2) (0) | — k,, F2r-2#-2) (zr) | 


C{F(x) } -f F(x) cos kyxdx = f.(Rn), 2 =1,2,---, (39) 
0 
where k, are the roots of tank r=(h/k), k>0, h¥0 
1§ Poe T \ a a : P 
{fe(Rn) } = >> N(Rn)fc(Rn) cos k,x in (0,7) and 
1 


1/N (Rn) = (w/2) + (4hq)~} sin 2har 
(Fen(x)} eae 


s=r—1 
4 > (- [cos kaw |FO r—28-1) (7) 4 fF Cr-2s-2) (77) | — Fr-2s 1) (0) | 


s=0 


C{F(x)} -f F(x) cos k,(m — x)dx = fc(Rn), n=1,2,--- (40) 


0 


a 
“- 
~~ 
3 
| 


, = z. N(Rn)f-(Rn) Cos k,(w — x), in (0, x) 
1 


k, and N(k,) as in (39) 


C{FO(x)} = (— ba) fe(bn) 
s=r—1 Pe 
." bs (—1)*+1, [cos kam [FO 28 1)(Q) — AF @r-2s 2) (0) | —F 2s (xr) | 


7 { F(x) } a F(x)(a sin k,x — k,b cos k,x)dx = f(kn), n=1,2,--- (41) 
i 
0 


where k, are the roots of tank t=k(bc—ad)/(ac+kbd), 
k> 0, ad — bc ¥ 0 
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T-"{ f(kn)} = >> N(Rn)f(Rn)(a sin ky — kyb cos knx) in (0,2) and 
1 


1/N(kn) = 3(wa? + wh,b? — ab)'+ [(k,b? — a2)/4k,] sin 2ek, + (ab/2) cos 2rkp 


2 s=r—] 
T{F@)(x)} = (— ka)f(kn) + DO (= k)* { ky [aF@r-2"-2 (0) + bFr- 22-10) | 


s=0 


— k,b(c sin Raw + Rnd cos ky)" [cF Or-2*-2)(e) + dF Or-2-D(m) J}, 
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NUMERICAL METHODS FOR FINDING CHARACTERISTIC ROOTS 
AND VECTORS OF MATRICES* 


BY 
W. M. KINCAID 
University of Michigan 


The present paper treats the problem of finding the characteristic roots and 
vectors of a matrix (having linear elementary divisors). The emphasis throughout is 
on methods of getting numerical results in practical cases rather than on theoretical 
questions. 

Symmetric matrices are taken up first, and methods are discussed for finding (a) 
the largest characteristic root and corresponding vector, and (b) the other roots and 
vectors. All these methods are variants of the iteration process. They are then ex- 
tended to general matrices, with particular reference to the case of complex roots. The 
paper closes with a brief discussion of the solution of algebraic equations by means of 
matrices. 

Among earlier work on this subject, we may mention that of Hotelling' and 
Aitken.? Indeed, much of the material in the present paper is taken from Aitken’s, 
through some modifications have been introduced, particularly with regard to the 
determination of roots other than the largest. 

A recent paper by Fry’ takes up matrices in connection with the solution of alge- 
braic equations. The present paper, particularly the last section, is thus in a measure 
supplementary to Fry’s. Wayland,‘ on the other hand, gives methods for reducing the 
problem of finding the roots of a given matrix to that of solving an algebraic equa- 
tion; this matter is touched upon Sec. II. 

Recent work along these lines has also been done by Morris® and Head.® 

At this point the author wishes to acknowledge his indebtedness to Prof. W. 
Feller for many helpful suggestions made during the preparation of the manuscript. 


I. BAsiIc RELATIONS AND DEFINITIONS 


1. Nature of the problem: definitions and notations. The problem of finding the 
characteristic roots and vectors of a matrix arises naturally in the solution of a sys- 


tem of linear equations of the type 


a,x) + a;.x “pb Saab a ai,x™ as Ax) \ 
Ao1x) + dgox 4 +--+ + donx™ = Aax@ | 
| (1) 
. \ 
Ani X™ +dnox + ee) + danx™ =Ax™ J 


* Received Dec. 27, 1946. This paper is a condensed version of a thesis submitted in partial fulfill- 
ment of the requirements for the degree of Doctor of Philosophy at Brown University, October, 1946. 

1H. Hotelling, Psychometrika 1, 27—35 (1936). 

2 A. C. Aitken, Proc. Roy. Soc. Edinburgh 47, 269-304 (1937). 

’T. C. Fry, Q. Appl. Math. 3, 89-105 (1945). 

4H. Wayland, Q. Appl. Math. 2, 277-306 (1944). 

5 J. Morris, Aircraft Engrg. 14, 108-110 (1942). 

6 J. Morris and J. W. Head, Aircraft Engrg. 14, 312-314, (1942). 
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In general, the system (1) has non-trivial solutions for m values Ay, Ae, - - +, An 
of the parameter A, which are called the characteristic values or roots of the system; 
these are the roots of the determinantal equation 


a4, —A ao} °° Ain | 
| 
a2 do2 — A+++ don 
| = 0. (2) 
| Gni An2 oo i r 
. . . . . . a ‘9 \ 
Corresponding to any characteristic value \; there exists a solution (x, xe. - - el) 
which may be termed a characteristic vector of the system. 
1 ( . . 
The vector (x{”, - - - , x{”) may be regarded as an n by 1 matrix. Such a matrix 
is called a column vector. Similarly, a matrix having but one row is called a row vector 
and will be written in the form jy, y®,---, yo} In general, if V is any given 


column vector, we shall denote by V’ the row vector having the same elements. 
If C= [ce] is an m by n matrix and D= [d;;| is an n by p matrix, their product CD 
is defined to be the m by p matrix E having elements of the form 


ei; = > cide; (¢=1,2,---,m;j=1,---, p) 
k=l 


Note that DC is not necessarily equal to CD. In particular, a row vector R and a 
column vector K (with the same number of elements) have two products, a scalar 
product RK and a matrix product KR. 

The system (1) may now be written in the form 


AX = XX, (3) 
where 

ay G12 °° °* Gin x) 
do, Go2°°** Aon - x?) 

A= X = ’ (4) 
Qn2  Gn2*** Ann xin) 

° . , (1) 
By a natural extension of meaning, the number \; and vectors X;=(x;", -- -, xi”) 


satisfying (3) are called the characteristic roots and characteristic column vectors of the 
matrix A. (Note that A is necessarily square.) Also any row vector Y such that 


YA = XY (5) 


for some value of \ is a characteristic row vector of A. Clearly the values of \ for which 
(5) is satisfied must be the same as before, for (2) must hold in both cases. 

In the following we shall suppose either that A has n distinct roots, or that to each 
root of multiplicity r there correspond r linearly independent characteristic vectors of 
each type (i.e., A has linear elementary divisors), as is almost always the case in 
practice.* The sets of characteristic roots and corresponding vectors will be denoted by 


* For a discussion of the case of non-linear elementary divisors, see the paper mentioned in Foot- 


note 2. 
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Mi, Ae, =) Ans X1, X2,°°*, Xnj Vi, Vo, - +--+, Vu. If any root A; is simple, the 
vectors X; and Y; are uniquely determined except for a constant factor; if two 
roots are equal, any linear combination of two corresponding characteristic vectors 
is a characteristic vector. In the following we shall understand by X,, Xi, Xi, 
(¢=1, 2,--+-+,m) the ith root and corresponding column and row vectors of some 
given matrix A of order n. 

Since (A; Y;)X;= Y;AX;= Y;(\iX;), we have (A;—A,;) Y;X;=0, and it follows that 
the characteristic column vector X; is orthogonal to the characteristic row vector 
Y; if \;#A;. We may clearly choose the vectors corresponding to any multiple 
roots in such a way that Y;X;=0 whenever 7#j, and we shall suppose this done in 
future. 

In what follow, unless stated otherwise, we shall assume that all the roots are 
distinct and shall take \, to be the root of largest absolute value. 

2. Some important relations. Since the X;’s are linearly independent, any column 
vector V of order m may be expressed in the form 


V = aX, + a2X2e + °° ° + anXn, (6) 
where the a,’s are constants. Multiplying (2) by A, we get 
A V = aA Xi + ios + haz = 43\1X1 + 2 -+- ra; ee (7) 


and by induction 
ue k k k ; 
A*V = Q4AiX1 + GoA2X 2 + al aa re ly a (8) 
for any positive integer k. More generally, if P(A) is any polynomial in the variable 
A, we have 


P(A)V = a,P(A1)X1 + a2P(d2)X2 + +++ + anP(An)Xn- (9) 


By setting V=X;, we see that the matrix P(A) has the characteristic vectors 
Xi, Xo, --+, Xn, and the characteristic roots P(Ai), P(A2), - - - , P(A,). (The same 


relations, of course, hold for row vectors.) 
A further useful relation may be derived from equation (2), viz., 


G11 + Goo + °° $+ Onn = Ar tA2t°°- +An, (10) 


i.e., the sum of the roots is equal to the sum of the elements of the principal diagonal. 

3. Symmetric matrices. A square matrix is symmetric if it is unchanged by turning 
its rows into columns and vice versa. The characteristic row vectors of a symmetric 
matrix have the same elements as the corresponding column vectors (Y;=X,’ 
for <=1, 2, - - -, ”); hence the characteristic vectors may be regarded as a single 
mutually orthogonal set. 

It follows that all the roots of a symmetric matrix are real, for two conjugate 
complex roots would correspond to two characteristic vectors having conjugate 
complex elements, and such vectors could not be orthogonal. 


II. THE METHOD OF ITERATION 


We now take up the problem of finding the characteristic roots and vectors of a 
given matrix A. For the sake of simplicity, we shall confine ourselves at first to the 
case where A is symmetric. 
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4. Expansion of the determinant. The most obvious procedure for solution would 
be to expand the determinant in Eq. (2) into a polynomial and then solve the result- 
ing algebraic equation for \. This method is sometimes the most expeditious, espe- 
cially if use is made of the techniques of expansion described by Wayland.‘ Never- 
theless it is necessary to consider other methods, for the following reasons. 

(a) Solving (2) gives only the roots; if the characteristic vectors are desired, they 
must be found by a separate process. 

(b) Finding the complex roots of an algebraic equation may be as difficult as 
finding those of a matrix directly. 

(c) In case approximate values of the roots are known, the work of expanding the 
determinant will not be affected, though that of solving the algebraic equation may 
be. On the other hand, we shall see in Sec. IV that such information can be used to 
great advantage when the problem is solved by iteration. 

(d) The methods to be discussed are better adapted to machine calculations. 

We therefore take up direct methods based on the prccess of iteration. 

5. Iteration. This method consists in multiplying a suitable vector V repeatedly 
by A. Representing V in the form (6), and considering equation (8), we see that if 
a;~0 (i.e., if V is not orthogonal to X) the first term on the right will predominate 


more and more as k increases. In fact 
AV = i(a.X1 + O(A2/1)*), (11) 


where Ag is the root of second highest absolute value. Similarly 


k+1 : ve 

AFHYy = A1 (a,X1, + O(Ae/A1) **?). 
Thus the ratio of corresponding elements of A*V and A**!V approaches \y as k>&. 
And if A*V is divided by one of its elements, the resulting vector tends to a multiple of 
X, (which, of course, is determined only up to a constant factor anyway). These 
points will be clarified by the following concrete case. 

Example 3. 

6 1 -1 519) 
1 4 O-—2 | 33 


~ 
un 
— 
on 


wn 
| 
bo 
wn 
~ 
oo 


(The extra column on the right is the sum of the others; its purpose will be explained 
in the following.) 
Choosing 


V = (4, —1, 2, 4), 
(more will be said later about the method of choice), we obtain 


AV = (33, —8, 18, 32), 


A*V = (268, —63, 145, 269), 
A’V = (2207, —522, 1222, 2193), 
A‘V = (18077, —4267, 9980, 18161). 
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Successive ratios of the leading elements are 
8.25, 8.12, 8.235, 8.191. 
(The true value of the root is 8.22557331.) The vectors obtained by dividing the fore- | 
going by their leading elements are (with the true vector for comparison) 
(1.000, —.242, .545, .970), 
(1.000, —.235, .541, 1.004), 
(1.000, —.237, .554, .994), 
(1.000, —.235, .552, 1.005); 
X, = (1.00000000, —.24073464, .55955487, 1.00862094). 


| 
bo 


| 
bo 


6. Use of powers of the matrix. It will be noticed that the convergence is not very 
rapid, and the multiplication might have to be repeated many times if A, or the 
components of X, were wanted to many figures. The process may be shortened by 
first squaring A repeatedly: 











P47 4 8 1777 76) 
: 4 21-11 -97 5 
Y=) og -41 27 12 | 361’ 
117 -9 12 42.) 62 
T2578 31 752 15737] 4934 
‘ee 31 659 —604 —631 | —545 
752 —604 1058 1063 | 2269 | 
1.1573 —631 1063 2278_} 4283 | 
A’ = +: 


(The powers of a symmetric matrix are of course symmetric also.) 

The sum columns are useful as checks at this point, for the sum column of A? 
should be the same as the result of operating with A on its own sum column, and 
similarly for the higher powers. That this is so can be seen by noting that the sum col- 
umn of A is equal to AH, where H=(1, 1, 1, 1). Then A*H =A(AH), AtH=A*(A7A1), 
Gt. 

We can now multiply V directly by some higher power of A and thus save many 
stepts in the iteration. For example 

A*(A4V) = A8V = (82542442, —19739077, 46035155, 83107096) 
A(A8V) = A®*V = (678801708, — 162628058, 379028193, 683495447) 


(5581640338, — 1338701418, 311703720, 5623793099) 


A(A®°V) = A*(A®8V) = AV 
These are proportional to 


(1.000, —.2391, —.5577, 1.0068), (1.000, —.2396, .5584, 1.0069), 





(1.0000, —.2398, .5586, 1.0076), 
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and the ratios of the leading elements are 
8.224, 8.224. 


Note, however, that forming each power of A requires about as much computation 
as m vectors A*V, so this procedure should not be carried too far. 

For example, to compute AV by operating directly with A on V sixteen times re- 
quires 16” multiplications; forming A? first and then operating eight times on V 
requires 3n(m+1)+8n multiplications. Similarly, using A‘ requires n(n+1)+4n; 
using A®, $n(m+1)+2n; and using A, 2n(n+1)+mn. For n=4, it will be seen that 
working with A‘ is the most economical. 

7. Finding the smaller roots. So far we have obtained only \;. What about the 
other roots? 

Clearly these can be found if we can determine the products AyAz, AideAs, etc. 
We shall now see how to modify the iteration process in order to do this. 

Suppose we have two sequences having kth terms of the form 


Bu = bir + beds +--+ + bade 
(k= @, 1, 2,--+). (12) 


, 


e & rk 7k 
Bi = bi\1 + Bore + lle + drrn 


Then we define 


| Be Brya | 
Bo. = r | 
By Brat 
= Apheldsts ~ BO, — Xe) + Ae — BIO — OD >, (13) 
whence 


Bo,441/B2,% = AiA2 + O(As/A2) * (k =0,1,2 :-- ). (14) 
To be concrete, if we have a sequence of vectors of the type V, AV, A?V,--:- 
A*V,--+we could write B,.=(A*V)™, By=(A*V)®, where (A*V)® is the ith 
element of A*V(¢=1, - - - , m), so that 
(ARV) (ARH) 
—_ (A kV) @) (A k+1//) (2) 
(of course, any other pair of elements would work equally well). Consider the follow- 


ing case. 


Example 2. 


1 1 4 -2 


‘It is found that \;= —8.075320860. Operating with A upon the vector V below 


yields 
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V = (5, —9 —2, 8), 

iV = (—12, 5, 75, —38), 

12V = (188, —258, —417, 393) 

18V = (—225, 1383, 3970, —2896), 

14V = (10836, —13423, —30093, 24280), 

15V = (—84466, 101136, 249340, — 193191). 
Then 

ed Oe Me Cedi Oe 
8 —38 —38 393 


Values of successive determinants and their ratios are 
— 524 13667 — 355474 9242272 — 240278437 


— 26.08 — 26.01 — 25.9985 — 25.9978 


Thus AjAz = — 26.00 and A, =3.219. 
On the other hand, we may take 6j=(;4:, in which case 


| Bi Br+t > 
Bon = (15) 


| Bi +1 Bro | 
That is, we can set up a sequence converging to \jA2 given any single sequence of 
the form (12) 
If we wish to find ); ,we select three sequences of the form (12), set up the cor- 
responding three-rowed determinant 63, and get as before 
B3,n+1/Bs,% = ArA2d3 + O(Ag/Az) *. (16) 


Note however, that the characteristic vectors Xo, X3,---,X, cannot be found 


by this means. 
III. DEvICEs FOR IMPROVING CONVERGENCE 

If we employ the method of iteration in the crude form, described above, we find 
that the convergence often leaves much to be desired in terms of time consumed. In 
this section we consider various ways of improving the efficiency of the method. 

8. Scalar products. Consider the scalar product (see Sec. I, 1) (V’A*)-(A'V), 
(k, =0, 1, 2,3,---). 
This may be written (in view of (6) and the orthogonality of X,’s) 


l 


V'AVIV = V'AMY = ap, XX, + ade XiXe+ +++ +ade XX (17) 


In particular 


V'AF-1.AFV = ashy X1X1 + aed 


ce 


sey UES 2k— S 
X 2X 2 + iain + anAn X nX my 


2k- 
2 


and 
2k__» i >. a 
= ay; X1X1+ aero X oXe+ > * + Anda XnXn, 


— 
— 
— 
> 
os 
‘ 


so that 
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V’A*-A*V 


V/AE. ARV = Ai + O(A2/A1)?*, (k = 1, 2,3,---). (18) 


It follows from (18) that by taking suitable scalar products the number of steps 
required to find \; to a given degree of accuracy is approximately halved. 

To illustrate this method, let us apply it to the vector sequence of Example 1. 
Working with the last two vectors, we get 


V’A'8V = 1098048041, 
V’A¥V/V'A8V = 8,.2255708, 
V’A*V = 9032871992, 
V’A®V/V’AV = 82255717. 
V’'A®V = 74293955655, 


(The last figures of the scalar products have been omitted.) The improvement in the 
approximation is evident. 

The method can also be applied to finding the products of the roots, by making 
proper use of (15). 

9. Choice of initial vector. So far nothing has been said about the choice of the 
initial vector V. Clearly, the best V is one that is already “near” to Xj, i.e., one having 
a, (equation (6)) large in comparison with the other a’s. A few suggestions may be 
helpful. The sum columns of A and even more of A?, A, etc., may be good indications, 
though sometimes they are misleading, as when H happens to be nearly orthogonal 
to X,; in this case we may be led astray to one of the other vectors. The other 
columns of A, especially those having the largest elements, are also good first guesses. 
In any case, the first V’s should have small integers for elements and should be tested 
by trial multiplications by A. 

10. Reduction of leading element to unity. Once V has been selected and the 
iteration is under way, it is clear (since only the ratios of the vector elements are sig- 
nificant) that the approximation to X; is not affected if each vector A*V is divided 
through by one of its elements—say the first if this is not zero. Then it will be clear 
from the extent of agreement of successive vectors how far the approximation has 
progressed, and time can be saved by retaining in the elements only as many digits 
as appear to be correct. (The choice of V is actually slightly altered at each step, but 
this is of no importance.) The following example shows the procedure. 


Example 3. 








{= 0-3 -5 i1 1 
3 5 i -3 —-1 
—1 0 1 —!] <i 


Casual inspection suggests the vector 
V = (1, 1, 1, —1, 0) 


as a first approximation of X,, The vector AV resembles V sufficiently to indicate 
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that we are on the right track. Operating repeatedly with A and dividing through 
by the leading element each time yields the sequence below: 
AV =(—7 —10, —9, 12, 1)(69), 
(1.00, 1.46, 1.28, —1.70, —.22)(—9.80), 
(1.000, 1.463, 1.296, —1.724, —.224)(—9.874), 
(1.600, 1.469, 1.298, —1.722, —.227)(—9.877), 
(1.000, 1.468, 1.301, —1.725, —.227)(—9.884); 


X, = (1.000000, 1.469801, 1.302061, —1.724997, —.228106). 


The figure in parentheses after each vector is the number by which the elements 
of the next vector have to be divided to reduce its leading element to unity. These 
numbers converge, of course, to Ax. Note that while computing errors may slow up 
the convergence, they will not invalidate the final result. 

To obtain a more accurate value of \i1, we should operate with A on the last 
vector in the sequence and take scalar products as before. The number of correct 
significant figures thus obtained should be about double that of the last number in 
parentheses. 

11. Aitken’s 6? process. Since the rapidity of convergence in all the above 
processes depends essentially on the ratio | )z, A: | , where dz is the root of second larg- 
est absolute value, the convergence is slow if th’s ratio is near 1. In such cases the 
convergence may be speeded up by using the 6? process, due to Aitken, which we 
now explain. 

Our approximating sequences are all of the type 


k k Fr 
de = Ar + C191 + C292 + °° (k= 0,1,2,---;1>]q1| =2lqe|2---) (19) 
where the g’s are of the form \2/Ai, As/A1, (A2/A1)?, etc. 
Writing 
| de Pk+ 
Ses = | : 3) (k =0,1,2,---), (20) 
Piri Pe+2 | 


we readily see that 
Yor = Ala(l — q)°41 + (1 — 92)°q2 fees] €1¢29192(41 — q2)?+--- 
On the other hand 
8 bisr = dete — 2pesi + oe = al — q1)°q1 + ¢:(1 — q2)"q2 dil 

Defining 

G2,k = Y2,4/O'br+1 (k = 0,1,2,---), (21) 
we see that 

$2.4 = 1 + O(q). (22) 


Thus, if | go| <| ql, the ¢2,’s will form a sequence converging to di more rapidly 
than the sequence of ¢x’s. Applying the same process to the ¢2,,’s will then result in a 
still more rapidly convergent sequence, and so on. 
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However, the remainder term in the sequence obtained by applying the 6* process 
twice is of the form O[(q3/q:)*|+O(qi). Thus the improvement in convergence will 
be slight if |gqi| and | g2| are nearly equal, especially if g.~—qu. 

12. Extension of the 6* process. For this case, we shall develop a modified pro- 
cedure which proves more efficacious. 

We define a second derived sequence 


3,k = W3,%/X3,k (k= 0,1,2,---), (23) 


/ $00 (24) 


where 


| py din. Pk+e | . 
| Wo,nqa Yo,k+2 
Var =| brit Orie Geis | = 
| | Werze Wenzs 
| Pr+2 Peis Pez | 





and 


| B*pisr SH hrre | 
a ; - (25) 


5*bi42 O-br+3 
It can be verified that 
k 
o3,% = A1 + O(9s). (26) 

It should be noted that if the elements ¢; of the original sequence are replaced by 
¢@. —c, where c is a constant, all the sequences will converge to \;—c. Thus if the first 
few digits of the ¢,’s are the same for all k, we may simply ignore these digits in 
computing the derived sequences. 

Example 4. 


2 1-1 3 1 | 
| 
1 0 0 -2 —1| 
A= V = 
-—1 0-3 5 3 | 
3 —-2 «5 —2 —3 | 


Carrying out the foregoing procedure, we form the following table: 


Arty dx 5px V2, 2,41 3,k-2 
20 — 8.550000000 — 8.59 
—171 — 8.602339181 .066643012 — .006071541 1105433 — 8.59 
1471 — 8.588035350 — .018447694 .001683322 1248386 120357 
— 12633 — 8.592179213 .005452930 —.000497222 1184409 120331 
108545 — 8.590870146 — .001769071 .000161358 1210537 


eo °) 


.591330150 


— 932496 _ 
8011381 
(In the above, ¢,—8.5, rather than ¢,, is used in computing the later columns.) 
Further derived sequences, showing still better convergence, could be defined in 
a similar fashion. However, the extra computation required and the rapid accumula- 
tion of rounding-off errors make such a procedure of doubtful value. 
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The same process can be used in finding the elements of X, (see the paper quoted 
in Footnote 2). 

IV. FINDING THE REMAINING ROOTs AND VECTORS 

In Secs. II and III we saw how all the roots of A could be obtained, though with 
limited accuracy for the smaller roots in practice. Now we consider methods of 
getting both roots and vectors at once, at the same time improving the accuracy of 
the former. 

13. \,-differencing. Perhaps the simplest method is that of \,-differencing the 
sequence (8). That is, we form the differences 


k k 
p9PF A w ofa = tdi + td = DS ++ HD — WA =D 


If X, is larger numerically than the other roots and a2;#0 the first term predominates, 
and we find X2 and dz, much as we originally found X; and ),. 

The disadvantage of this method is that nearly equal numbers have to be sub- 
tracted at each stage with a resulting loss of significant figures. The better the con- 
vergence of the original sequence is, the worse this becomes. Needless to say, the 
difficulty is increased if we try to find X3, X4, etc., by repeating the procedure. 

14. Deflation.!* We may also proceed by deflation, which consists in replacing 
the matrix A by a matrix A; having the same characteristic roots and vectors as A, 
except that A; is replaced by zero. Applying the method of iteration to A; will then 
yield \, and X2. Deflating again give \3 and X3, and so on. 

We set up the matrix A; as follows: 

Ai 


A— XX. (28) 
X{X, 


Note that X/ X, is the scalar product (sum of squares) while XX is the matrix 


product 


— 





(1)2 (1) (2) @ 4 
x1 X1 - 41 %1 
(2) (1) (2)2 (2) (n) 
xX, 44 x1 - 4, 
(n) (1) (n) (2) (n)2 
= aes %1 %1 x1 a 








That A, has the desired properties can be seen by substitution: 


At 
i XX X; = 


— 1.50364702_j— .08483344 


1,X1 = AX, sai AX) — AX] = QO, 
Xi xX, 
r r Ai r r . 
1X; = AX "a = gee Rane Xi = MX; (i - 1). 
Xi X, 
Let us apply this procedure to Example 1: 
XX, = 2.38837102, Ai/X1 Xi = 3.44400985 
2.55599015 1.82909247 —2.92711248 — .47370045 . 98426969 | 
; 1.82909247 3.71445974 .46392273 —1.16375997 | 4.84371497 | 
; —2.92711248 .46392273 — .07832518 3.05627400 .51475907 | 


— .47370045 —1.16375997 3.05627400 
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Starting with the vector V; below and proceeding as in Sec. III, 10, we get in succes- 
sion the vectors 


Vi = (2, 2, —1, —1)(12.216), 
(1.00, 0.96, —0.65, —0.40)(6.4040), 
(1.00, 0.87, —0.57, —0.46)(6.0337), 
(1.00, 0.88, —0.64, —0.42)(6.2379), 
(1.00, 0.85, —0.60, —0.45)(6.0802) 


whereas the vector and root to four places are 
X-2 = (1.0000, 0.8681, —0.6198, —0.4404), A» = 0.16666. 


We see that the difficulty of losing significant figures by subtraction is largely 
avoided, while the convergence toward X; is fairly rapid. On the other hand, the 
extra figures in the elements of A; make the further computations somewhat labori- 
ous, and in any case the accuracy attainable in X2 is limited by that of \; and Xi. 

Some saving of time may be affected by noting that 
Ns 


. 
A,=A*f— commana XiXi1, 
xi dz: 
whence 
k 

i ew 
—— XXi J, 
XX, 


i.e., the vector ASV can be obtained from A*V without setting up the matrix Ai. 
But this again brings a loss of significant figures; thus this device is of limited ap- 


SY = AV = 


plicability. 

15. Use of matrix polynomials. In view of what has been said above, it would 
seem that there is room for a different approach. Now, upon reviewing our work up 
to this point, it will be seen that we have tacitly assumed that we had no prior notion 
of the location of the roots of our matrix A, and we have indicated no way of using 
such information if we had it. However, in a practical case growing out of, say, a 
vibration problem, we could usually make at least a guess at where the roots lie. In 
any event, we could get an approximation to the largest root by a couple of steps 
of iteration, and we shall see shortly how the other roots might be estimated. 


Suppose, then, that we have approximate values Ay, 7, ---,An of the m roots 
Ai, As, - - » » An» Now consider the polynomial. 
P(A) = (A — AZ )A — Az) +++ A Ad) (28) 


and the corresponding matrix 
P(A) = (A — dg)(A — AZ) +++ (A — Ae). (29) 


It follows from equation (9) that P(A) has the characteristic roots Pi(\1), 
Pi(X2), - > >, Pi(\n). But by definition 


Pi?) = Pads) = +++ = Pids) = 0, 
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and since P(A) is continuous one would expect Pi(Az), Pi(As), - ~~, Pi(An) to be 
small, how small depending on how good the approximations were. On the other hand, 
the roots being supposed distinct, Pi(Ai1) would not in general be near zero, and 


therefore 
|PiQs)| >|] PiQ)|, i 1. 
If we now choose a vector V of the form (6), with a,;#0, we have by (9) 


P,(A)- V = a,P3(A1)X1 + d2P;(A2)X 2 aa * 3 + AnP (An) Xn- (30) 


The first term on the right will predominate strongly, so that if we proceed by itera- 
tion, i.e., multiply repeatedly by P(A), the convergence to X, will be very rapid. 
(Of course it is all the better if V is already an approximation to X.) 

Once X; is obtained, \, can be determined to the same degree of accuracy by 
multiplying X, by A, and to greater accuracy by forming scalar products. 

Note that nothing is really altered if \; is replaced by one of the other roots. For 
example, to find Xz: and X2 we could use 


Pia mm 14 AI ~ GA <)> = 2} 
Example 5. Suppose that for a given fourth-order matrix we have 


Ai ~ 6.1 = hj, ho ~ — 4.3 = do, A3 ~ —2.6=)N3, he~ 0.4 = rf. 


Then 
P,(d) = (A + 4.3)(A — 2.6)(A — 0.4) = A3 + 1.342 — 11.94 + 4.5, 


where the coefficients have been rounded off to one decimal place, as there is clearly 
no advantage in retaining more. In fact, we should in practice make a further simpli- 
fication, and replace P;(A) as written above by a polynomial like \*+A?—12A+4 
or 2\°+ 3A? — 24+-9, whose coefficients are small integers. In this way, the subsequent 
computations would be simplified without slowing down the convergence much. 
Where to strike the balance in practice is largely a matter of experience and “feel” 
on the part of the computer. 

16. Application of orthogonality. Let us return to Sec. I, 2, but suppose that IV, 
instead of being a completely arbitrary vector, is orthogonal to X,, i.e., ai=0. (We 
shall see later how to select such a V.) Then 


iF 6 ei + tbs & + +> +O, 
P(A)V = d2P(r2) V2 + a3P(A3)X3 + dell +- GF th.) Aw 


(31) 


To see the application of this, consider Example 5 again and suppose that X; 
has been found to the desired accuracy. From what has been said before, our first step 
in gett X_ would be to set up the polynomial 


P2(d) = (A — 6.1)(A — 2.6)(A — 0.4). 
However, if we are to work with a vector V orthogonal to Xi, we might just as well 
use the simpler polynomial 


Q2(A) = (A — 2.6)(A — 0.4). 


aw 
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The fact that Q2(A1) is not small compared with Q2(Az) has no effect on the iteration, 
for 


Q2(A)-V = aeQo(A2)X2 + a3Q2(As)X3 + adQo(Aa) Xa. 
- a ; vi”) 


How can we select such a vector V? Consider any vector Vo=(v’, 0%; 
Then the vector , 
(i) G) (2) (2) (n—1) (n—1) 
( (1) (2) (n—1)  %1 +% 4% +°°* + *1 ) 
% 3 — —_______—_—_— 


YM »°** 4% — —— ae (32) 
“ 
will be orthogonal to X; and can serve as V (provided x!” 40). 

Note that we have altered only the last element v” of Vo. We could, of course, 
have altered any other element v instead, provided x{? #0. In general it is best to 
alter an element corresponding to one of the larger elements of X; so that V may 
differ as little as possible from Vp. 

In order to avoid the accumulation of rounding-off errors, it is generally desirable 
to repeat the orthogonalization at each step of the iteration. Also, it is well to de- 
termine the elements of X; to one more significant figure than is required in those 
of the other vectors. (Of course we could use any of the vectors X; as a starting point.) 

17. Obtaining the first approximations. It is well to give some attention to the 
question of getting rough values of the roots when these are not available before- 
hand. Since the procedure is rough and ready rather than fixed, it is best shown by 
an example. For this purpose we return to the matrix A of Example 3. 

We have already seen that 


Ai = — 9.886, X; = (1.000, 1.470, 1.302, —1.725, —0.228). 


Looking back, we see that one or two steps in the iteration would have sufficed to 
give us one-place accuracy, which is all we need at this stage. 
Referring to equation (10), we see that 
Ar tAs+tAsHtAG+As = —24+0-—-5-3-1=-IH11, (33) 


so that A must certainly have a root dz greater than —1. It follows that the matrix 
A, =A+6 will have a root \2.+6>5, which is greater in absolute value than its root 
i +6. Trying out simple vectors like (1, 0, 0, 0, 0), (0, 1, 0, 0, 0), etc., we see that the 


combination 


V; = (0, 1, 0, 1, 0) 


looks hopeful. In fact 
Aw, = (1, 11, -2,8, —1), Ai = (7, 110, —28, 81, —16), 


and by taking scalar products we get \2+6~10.2, \2~4.2, while A?V, will stand as 
a first approximation to X2. 

To proceed further, we note that \i+3~—6.9, \2+3~7.2. Therefore the matrix 
(A +3)2—50 will have two roots near 0 and three between —50 and 0. This matrix 
could be used to obtain another root A; of A. However, it probably has two or three 
roots near —50, which would be difficult to separate, so we used instead the matrix 


Ay=(A2?+6A +9) —35=A(A +6) —26. 
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After a few trials we decide on the vector V2=(2, 0, 0, 1, 4). Then A2V2=(—70, 
—4, —6, —27, —128), and taking scalar products yields 


(A3 + 3)? — 35 ~ — 32.5, (A3 + 3)? ~ 2.5, As +3 ~ + 1.6, 
so that \x~ —1.4 or —4.6. Since A V2=(—5, 1,5, —1, —7) and Ai V2=(7, 1, 5, 5, 17), 


the value \3~—1.4 must be the correct one. 
To get approximations to \y and \;, we take the sum of the principal diagonal 


terms in A>: 


: 2 
> A: + 6A; — 26) = — 55, 
i=1 


so that, from (33) 
Mi he + de the t+ As = 141. (34) 


Combining (33) and (34) with the known values of Ai, Ae, and A3, we get Ay~0.9, 


18. Examples. Two more examples show how to carry through the process. 
Example 6. Find to six decimal places the roots of 


2 2 0 + 8 


6 «=f =f 94 3 
ee wt @ oeTuy 
‘$s 2 ef 8 


Using the methods of the last example, we conclude that \;~7.06, \»x~—4.19, 
As~0.45, Ay~w—2.32, Xi~(106, 67, —38, 100), Xe~(53, 53, —41, —106)= V2. 
We next compute X, to greater accuracy in order to obtain ), to six places. Since 
we want \, to seven significant figures, four figures in the elements of X,; should be 
sufficient, as the formation of scalar products about doubles the number of signifi- 
cant figures, as explained above. In order to avoid rounding-off errors, both here 
and in the later calculations, we carry an extra figure. 

If we were to follow rigidly the procedure outline above (paragraph 18), we should 
begin by writing the polynomial (A+4.19)(A—0.45)(A+2.32). However, we note 
that the simple /inear matrix polynomial A,=A,+2 has very favorable root ratios, 
and since the number of figures required is not great anyway, it seems worthwhile to 
go ahead with this matrix rather than taking the time to set up a more complicated 
expression. Applying A: to our first approximation and dividing by the leading 
element at each step yields successively 

(106, 67, —38, 100)(958), 

(1.000, .644, —.358, .941)(9.052), 

(1.0000, .6435, —.3582, 19423)(9.0562), 
(1.00000, .64360, —.35826, .96206)(9.05544), 
(1.00000, .64360, —.35826, .94213)(9.05572), 


(9.05572, 5.82825, — 3.24438, 8.53158). 
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The last vector was obtained by applying Az once more without division. Taking 
scalar products of the last two vectors, we get 
2.4301801245 
22.0069707442 9.055695307 
199.2884216617 9.055695307 
We conclude that A; = 7.055695307, to greater accuracy than was required. 
We could now find ), in the same manner by starting with A —2, say. However 


the convergence would be rather slow because the ratios of the roots would be less 
favorable than before. Instead, we proceed as in Example 5. Since 


(A + 2.32)(A — .45) = A? + 1.87A — 1.04, 


we set up the matrix A;=A?+2A—1. Now if we start with a vector orthogonal to 
X;, we shall get a sequence converging to X2. We therefore start with the vector V2, 
but modify its last element so as to make it-orthogonal to X:. The modified last 


element is 
— [53(1.00000) + 53(.64360) — 41(— .35826)]/.94213 = — 108.05. (35) 
We now apply A; to the modified vector several times, repeating the orthogonal- 
ization process at each step. We get successively 
(53, 53, —41, —108.05)(417.90), 
(1.000, 1.130, —.739, — 2.11438) (8.21364), 
(1.0000, 1.1257, —.7384, —2.11122)(8.19976), 
(1.00000, 1.12581, —.73843, —2.11130)(8.20028), 
(1.00000, 1.12579, —.73842, —2.11128). 
Applying A,=A —2 (of which \.—2 is the numerically largest root) and taking scalar 


products yields \,.—2 = —6.1937207, \»x= —4.1937207. It follows by direct calcula- 


tion that \3= 0.464791, \4y= — 2.326766. 
Example 6. Find to five significant figures the characteristic vectors of the matrix 


> 2 + 

1-1 -1 3{ 2 
B= | 

0-1 8 —2| -3 | 


First we compute B? and B*. 

Proceedings as in the previous example, we shortly conclude that X:~(76, 44, 
—30, 76), X2~(67, 84, —37, —129), \1~6.56, Ax~ —4.53, As~0.69, Aw — 1.72. 

Next we find X, to the desired accuracy. Since more accuracy is required than 
in the previous example, it is worth while to go to more trouble to secure rapid 


convergence. We have 


(\ — 0.69)(A + 1.72)(A + 4.53) = A? + 5.56d? + 3.48A — 5.39. 
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Thus B,=2B*+11B+7B—11 will have three roots near zero. We get successively 


(76, 44, —30, 76)(81506), 
(1.0000000, .5789512, —.3916276, .9944053)(1069.511549), 
(1.0000000, .5789483, —.3916329, .9944266)(1069.520900), 


(1.0000000, .5789483, —.3916329, .9944265). 


The rapidity of convergence is to be noted. (As before, the extra digits are retained 
to avoid rounding-off errors both here and in the orthogonalization that follows.) 

To obtain the vector X2, we write (A—0.69)(A+1.72) =A?+1.03A—1.19, and 
set up the matrix B,=B?+B-—1. As in the previous example, we must operate on 
vectors orthogonal to X;. Orthogonalizing at each step we get in succession 


(67, 84, —37, —129)(1048), 

(1.0000, 1.1698, —.5821, —1.915902)(14.877570), 

(1.0000, 1.1894, —.5893, —1.930148)(15.003080), 
(1.000000, 1.189683, —.589403, — 1.9303534)(15.0048880), 
(1.000000, 1.189687, —.589404, —1.9303561). 

The convergence could have been improved by setting up a cubic matrix poly- 
nomial as was done with X,, but the extra time needed for doing so would have 
more than offset that saved later. 

We could now obtain X; by starting with B+2, say, and working with vectors 
orthogonal to both X; and X2. The double orthogonalization, however, would be 
time-consuming, so we start again with (A+1.72)(A+4.53) =A?+6.25A+7.79. 

Use of B?+6B+8 suggests itself. However, a quick estimate of the roots indi- 
cates that this will not lead to very rapid convergence, and as yet we have not even 
a first approximation to X3. So we use instead B;=4B?+25B+31. 

Starting with a vector that turns out to be a rather poor guess and orthogonalizing 
with respect to X; at each step yields the following rapidly convergent sequence: 

(1, 0, 0, —1), 

(21, —10, 26, —5)(1036), 

(1.0000, —.5097, 1.2838, —.203265)(50.266140), 
(1.000000, —.516826, 1.293370, —.195347)(50.513110), 


(1.000000, —.516805, 1.293339, —.1953714). 


The vector X4 can now be determined from the fact that it is orthogonal to 
Xi, Xe, and X3. We set 
Xe = V+ 1X1 + 42X22 + 3X3, (36) 
where V is any vector not orthogonal to X/. Multiplying (36) by XY gives 
0=Xi/ V+a,X/ Xi, whence 
a, = — X{V/XiX, (37) 
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and similarly for a; and a3. Choosing V=(1, 0, 0, 0), we obtain a,= —0.4036424, 
a2= —0.1541062, a3= —0.3357978. Substituting in (36) and dividing through by 
the leading element of the resulting vector yields X,=(1.000000, —2.287233, 
— 1.741508, —0.359851). 

The characteristic roots can be found, if desired, simply by computing the lead- 
ing elements of BX,, BX2, BX3, and BX4. Greater accuracy can be achieved by form- 
ing scalar products as in previous examples. 

19. Summary. We now summarize our procedure. 

(a) Obtain rough values of all the roots. This may be done by operating on 
properly chosen vectors with matrices of the forms A-+n, (A-+m)*+n, etc. After 
all but two of the roots have been located in this way, the last two can be found by 
solving a quadratic equation. 

The examples indicate that three significant figures in the roots are amply suffi- 
cient at this stage, and one could often get along with fewer. Certainly a dispropor- 
tionate amount of time should not be spent in getting these first approximations, as 
their sole purpose is the saving of time in the later steps. Fortunately the question is 
beside the point in many practical cases, where rough values of the roots are known 
to begin with; this step can then be omitted. 

Naturally, if only a few figures are wanted in the final results, it may prove most 
economical to carry the work to a conclusion at this stage, rather than taking the 
time to set up the matrix polynomials of the succeeding steps. 

(b) Select a polynomial P(A) such that, for some ),, | PA,)| >| PA,)| , j/ #1, and 
operate on a suitably chosen vector with P(A) until X; is obtained to the desired 
degree of accuracy. 

P(A) will ordinarily be chosen as in the examples. It will be noted that one or 
two figures in the coefficients are enough for rapid convergence, and more should 
not be used as time will be wasted in computing the elements of P(A). It should be 
borne in mind that the coefficients of the higher powers are the most important. 

Care should be taken to avoid errors at this stage, as much time can be lost if 
P(A) is incorrectly chosen or written. 

An extra place should be retained in the elements of X; to avoid rounding-off 
errors, and it may be well to keep two such places if vectors orthogonal to X; are to 
be constructed in step (c). 

(c) Proceed similarly to obtain the other characteristic vectors, making use of 
the orthogonality of the vectors as in the examples given. If only the roots are de- 
sired, two of the vectors need not be computed. If all the vectors are wanted, the last 
one may be obtained from the others by the condition of orthogonality. 

(d) Once the vectors are obtained, the corresponding roots may be found to about 
twice as many significant figures as the vectors by the scalar product method. In case 
the order of the matrix is greater than four, it may be necessary to use a matrix of 
the form (A-+m)?+n to do this, getting first the square of the root and then finding 
the root itself by taking a square root. The last two roots may be obtained from the 
others by solving a quadratic as before. 

In case A has multiple roots, this fact should become evident while carrying out 
step (a). We can then proceed much as before, except that we should not, of course, 
expect to find unique vectors corresponding to the multiple roots. 
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V. GENERAL MATRICES; COMPLEX ROOTS 


We now approach the problem of the general (unsymmetric) matrix. Unlike 
the symmetric case, the characteristic roots are no longer necessarily real, and there 
are two sets of characteristic vectors instead of one, having less convenient ortho- 


gonality properties. 

Most of the methods of the preceding section remain applicable. However, we 
cano no longer make such effective use of the orthogonality relation. Also, since the 
elements of the row vector V’A need not be the same as those of the column vector 
A V, the method of scalar products no longer effects a great saving of time. On the 
other hand, new methods are required to handle complex roots. 

20. Real roots. The case where all the roots are real presents few difficulties, and 
the procedure is very similar to that for the symmetric case. Some use may even be 
made of the orthogonality properties. 

The same can be said if only two of the roots are complex, and the real roots can 
be found first, since the complex pair can then be obtained by solving a quadratic 
equation. The following example illustrates these points. 

Example 7. 


1-3 0 —2]-4] 
2 
6 | 


Find all the roots correct to six places, given that \y~2.49. A2~ — 1.77. 
We compute A? and A# and deduce that A3,4~ — 0.36 + 3.282. Since 


(A + .36 — 3.282)(A + .36 + 3.282)(A + 1.77) = A® + 2.45A? + 12.2A + 19.3, 


we can get \; by using the matrix A,;=4A*+104+49A+77. After a few steps we 
obtain 
X, = (—.3637407, .0336918, 1.0000000, —.2743017), 
whence (multiplying by A once) \; = 2.4868715. 
To find As, we make use of the fact that the characteristic row vector }2 is 


orthogonal to X;. We have 


(A + .36 — 3.281)(A + .36 + 3.287) = A* + .72A + 10.9. 


Writing 


oe 199 <r {Ory As 
-63 9% -10 36 61 
—54 37 183 57 223 | 
oe wae! 269 ‘Wey aus | 


Az = 10A? + 7A + 109 = 


we get the following sequence of row vectors by multiplying on the right by Az and 
orthogonalizing with respect to X; at each step: 





1947] FINDING ROOTS AND VECTORS OF MATRICES 339 


{0, 1, 0, 0}, 
{1.00, —1.56, —.46, —0.57}, 
{ 1.0000, — 1.5795, 0.26538, —0.5526}, 
{ 1.000000, — 1.579690, 0.2653589, —0.552692}, 
{ 1.0000000, — 1.5796883, 0.2653587, —0.5526927}. 


Operating once with A yields \:= —1.7684837. It follows that A34= —0.3591939 
+ 3.2840604z. 

How to find the characteristic vectors corresponding to such complex roots will 
be explained below. 

21. Complex roots. When the complex roots cannot be approached from behind 
in this fashion, a new method of attack is necessary. Since such roots come in conju- 
gate pairs of equal modulus, applying the matrix repeatedly to an arbitrary vector 
will not result in a convergent vector sequence (in fact such behaviour is sometimes 
the best indication that the largest roots are complex). 

However, if \;=71e"1, Ae =ne7*1, |Aa| = | Ae| > |As| > ---+,wecan still find the 
product \,;A;=r; by the method used in Example 2; indeed the-convergence should 
be relatively better than when |u| > | Aa| , since the loss of accuracy due to subtract- 
ing almost equal numbers will not arise. Having r7, we can find )\,; and ), as follows. 

First by operating repeatedly on a vector V with A (or better, with some power of 
A) we obtain a vector W of the form 

W = aX. + a2X2+ A, (38) 
where A is as small as we please. It can now be seen that 
6 habits ni 0 in 12,3 ) (39) 
r, Cos 0) = ————_______- A j= 1,2,3,-+-, 9). 
1 1 (A wy + ( ) Gj ’ 
Thus A; and do are determined to any desired accuracy. 
Furthermore, 


AW + i(r?W — A*W)/2r, sin 0; 





2 2 
AA2(aiX1 + 2X2) — ariXi — ashaX 
= aiiX1 + oseX2 + §————__——_—————— + 04) _ (40) 





= 2a;A\1X1 +. 0(A) 
so that X, and X_ can be found once \, and dg are known. 

The problem can also be approached by operating with a matrix of the form 
A-+ni, whose roots will ordinarily have distinct moduli, but in practice this method 
will usually prove more laborious than the preceding one 

If all the roots are wanted, we begin by locating them roughly, as in the sym- 
metric case. Then matrix polynomials can be constructed to secure rapid convergence 
to more accurate values. These points will become clearer in the following example. 

Example 8. 

1 -2 0 —4] -5 
3 0 1 2 6 
—1 3 -1 1 2 


1 0 4 0 5 








340 W. M. KINCAID [Vol. V, No. 3 


Find the roots of A to seven decimal places and its characteristic column vectors to 
seven significant figures. 
As usual, we write down A? and A?: 


—j «fags «87 =S? | af ui 143-63 
4 —3 7 —11 —3 13 —23 13 —54 —15] —79 
A?= Pir 
10 —1 8 9 26 8 4 27 —34 5 
—3 10 -—4 0 3 31 -6 14 28 67 
Thus >°4.,’,=0, >-4., A? = —4, and the largest root must be complex. 


Starting with the vector V =(0, 0, 1, 0), we operate repeatedly with A*, rounding 
off to three or four figures at each step to avoid accumulating useless digits. After 
several such multiplications we come to the sequence of vectors given below: 

(—458, 752, —68, —432), (—2974, 3046, 1220, —3176), 
(1588, 898, 1293, —1823), (— 7060, 572, 8778, — 8756). 


As we proceed, we form two-rowed determinants from the first and second ele- 
ments of successive vectors, and likewise from the second and third elements, and 
then take the ratios of the values of successive determinants. The terms of tehse 
sequences corresponding to the above vectors are 


841380 2166396 5431544 


bo 
wn 
~~ 
bo 
wm 
_ 


1224568 2482918 7143048 
253 251 
(The sequences of ratios above show convergence to r{, aside from the location of the 
decimal point, which is immaterial at this stage.) 

The agreement of the last numbers suggest that if we now operate on the last 
vector with A and form ratios of determinants as before, we shall obtain rj to about 
three figures: 

(— 7060, 572, 8778, —8756), (26820, —29914, —8758, 28052), 
(— 25560, 127806, —79752, —8212), (—248324, —172856, 480518, — 344568), 
19585 266316 3615550 
13.60 13.58 
25758 350503 4762747 
13.61 13.59 


Thus 47~13.6. Using equation (39) with W as the second vector above yields 
7, COS 0; ~ — 2.27, r, sin 6; ~ 2.91, 


Ar~ — 2.27 + 2.911, Ar~™ — 2.27 — 2.912. 


It follows that A3.4=2.27 + 1.967. 
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In order to get the roots more exactly, we note that 
(A — 2.27 + 1.961)(A — 2.27 — 1.96%) = dA? — 4.54 + 8.99, 
Thus the matrix A;=2A?—9A +18 will have two roots near zero. Rounding off the 


elements of the last vector above, we multiply it by A; and obtain ratios of de- 
terminants formed from the first and second, and second and fourth elements: 


(— 2483, —1729, 4805, —3446)(— 243759, 188294, 153735, —27083), 
(—6120173, 18490944, —8355979, —4216169), 
(53044664, 46504148, — 111097501, 75339351), 
(56799532, —40188604, — 39094634, 63788114), 
— 888993313 — 335494216 — 126545934 — 477320484 
3774 37719260 37719148 
111704493 421315396 158916507 599419524 
3771 37719131 37719148 


We now multiply the last vector above by A and form ratios of determinants, 
getting 
(56799532, —40188604, — 39094634, 63788114), 
(— 117975716, 258880190, — 74482596, — 99579004), 
(— 237420080, — 627567752, 869519878, —415906100), 
9962994304 1355011183 
13.6004404 
1311417783 1783585959 
13.6004405 
Taking rj =13.6004405, we obtain, by using (39), A1,2= —2.26774878 + 2.908222131. 
Using equation (40) with W as the first vector above gives 
2a,X, = (— 117975716 + 1736316362, 258880190 + 1392346112, 
— 74482596 — 2409073411, — 00579004 + 2206593057) 
X,; = (1.00000000, — 0.63822188 — 1.05732751i, 
— 0.74982611 + 0.938445732, 1.13604812 — 0.198391772). 
The remaining roots must be A3,4=2.26774878 + 1.956428667. 
22. Finding the remaining vectors. The vector Xs; cannot be determined as 
X, was, since there is no ready-made vector sequence available. However, the matrix 


polynomial 


II 


(A + 2.26774878 — 2.90822213i)(A + 2.26774878 + 2.90822213%) 
-(A — 2.26774878 — 1.956428667) 

(A* + 2.26774878A? + 3.31507144A — 30.8423824) 

+ i(1.95624866A? + 8.87337741A + 26.6082916) 


Agog 
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has three roots so near zero that a single multiplication of an arbitrary vector by A? 
should give us X; to about eight figures. To do this it ts not necessary to compute the 
elements of the matrix A», for by (9) 


A2V = A®V + 2.26774878A°V + 3031507144AV — 30.8423824V 
+ i(1.95642866A°V + 8.87337741AV + 26.6082916V). 
Taking 
2 ee 
we see that AV, A?V, A*V are the sum solumns of A, A?, A*. Then 
AsV = (— 174.3244445 — 90.14645597, — 96.7552001 + 73.9792701i, 
39.7292288 + 95.22219167, 59.5362211 + 76.84446467) 
X3 ~ (1.00000000, .26477276 — .56129590i, — .40277933 
— .33795067i, — .44932357 — .2084592272). 
(The second vector above was obtained from the first by dividing by the leading ele- 
ment.) As a check, we perform the same operations with 
V; = (1, 0, 0, 0) 
(note that A Vi, A?V, A*V,; are simply the first columns of A, A*, A*) and get 
(1.00000000, .26477273 —.56129590i, —.40277934 —.337950661, 
— 44932357 —.208459221). 


VI. SOLUTION OF ALGEBRAIC EQUATIONS BY MATRIX METHODS 


23. The companion matrix. Before closing we consider briefly the solution of an 
algebraic equation 


J) =A" — pr — pre? —--- — p, = 0 (41) 
by the use of the “companion matrix” 
i” pe eeeee Pa Pa] 
Ly Os «+ 0 
Cmi0- is< 26 0 (42) 








10. O-+-1 0 


(Bernouilli’s method). It is easily seen that (42) has the same roots as (41) and that 


the characteristic column vector of (42) corresponding to the root \;(¢=1, 2,--- , m) 
is 
i 2 aid —3 —n 
Lek ctw OSA kh (43), 


Thus the roots of (41) are determined if either the roots or the column vectors of 
(42) can be found. 
24. Special methods of solution. While the methods employed in handling a gen- 
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eral matrix are applicable in this case, the special properties of C can be used to 
speed up the solution. First, equation (41) helps us in locating the roots, especially 
the real ones. Again, if \{ is an approximation of \,, the matrix polynomial 

v / 

A-—A 
(which may be obtained by synthetic division) will have only one large root, cor- 
responding to \;. We recall that to set up the corresponding polynomial P(A) in the 
general case approximations to all the roots are required. (However, (44) is not 
easily set up if \; is complex.) 

Moreover, suppose we have a vector V that is an approximation of X; and that 

D=P(C) is a matrix having P(A,) as its largest root. Then our next approximation, 
under our previous procedure, would be 


DV 


© (DVO ” 


W = (1, we, ws,+ ++ , Wn) 


Instead we may use the more easily computed 


e 2 n—1 
W* = (1, we, we, We ). 

This can be done because X; and W* are both of the form (43), and because w+, 
being obtained from. W, is an approximation to \;'. Unfortunately this process is 
not necessarily convergent (though in practice it generally is, especially when 

P(d,)| >| P()|, 742). However, we can always try it and then fall back on the 
standard process if necessary. 

Clearly, rather than determining W* from the first two elements of DV, we could 
determine it from the last two. It can be shown that this will lead to better con- 
vergence if \; is the smallest root of (41), while the first process will work better if 
\; 1s the largest root of (42). For roots in between we shauld use one process or the 
other depending on their relative magnitude. 

25. Examples. The following examples will make this procedure clear. 

Example 9. 


fA) ='-—- 4-3 =0. (46) 


Direct substitution shows that (46) has real roots between —1 and 0 and between 


1 and 2. 
If we take 2 as a first approximation to \,;, the matrix corresponding to (44) is 


C,=C*+2C7?+4C+4. 
The vector corresponding to \; can be written 
. ae 
Xi = (1, Ay ,A1 ,A1 }. 
The first element is taken as unity since ), is a relatively large root. Since \:~2, we 
start with the vector V=(1.00, .50, .25, .12). Operating with C, yields the vectors 
(1.000, .561, .315, .177), (1.0000, .5604, .3140, .1760). 


In forming these vectors, the first two elements of each vector were obtained in 
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the usual way, but the last two elements are simply the square and cube of the 
second element. Note that only the first two rows of C, are used in this procedure. 
From the last vector we get 


1 
Ay ~~ —— = 1.794. 
0.5604 
In order to get \; more exactly, we divide (46) by \—1.784 getting approximately 
Ae + 1.784? + 3.18 + 1.68. 


Thus the matrix C;=5C*+9C?+16C+8 has only one large root, corresponding to 
Ai (cf. (44)). Operating as before (using only the first two rows of C2) gives the 
rapidly convergent vector sequence 


(1.000, .560, .314, .176)(1.00000, .56043, .31408, .17602), 
(1.0000000, .5604256, .3410769, .1760167), 
(1.00000000, .56042566, .31407692, .17601677), 


whence A; = 1.7843580. 
Note that it was not necessary to approximate the other roots in order to find ),. 
The root Az can be obtained in much the same way. Starting with the matrix 


Cs=C?-C?+C-5 
and the vector 
V = (-—1, 1, —1, 1) 
we reach after a few steps the vector 


(—.334, .482, —.694, 1.000). 


(Since dz is the smallest root, we take the last element as unity here and use the last 
two rows of C; in the computation.) Thus \x~ —0.69. Dividing (36) by \+0.69 gives 
A? — 0.69A? + 0.48A — 4.33. 

Using the matrix 
Cy, = — 10C* + 7C? — SC + 43 


yields after a few more steps A, = —0.69250484. 
The complex roots are then A3,4= —0.5459266 + 1.45937791. 
Example 10. 


M—4+4=0. (37) 
Here the roots are all complex. The companion matrix is 


0 0 4 —4]0 
1 0 0 0} 1 
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By the methods employed in the general case, we get the rough values \y.2~—1.05 
+ 1.432, As 4~1.05 + 0.412. 
Noting that 

(A — 1.05 — .417)(A — 1.05 + .412)(A + 1.05 + 1.432) 

= (A* — 1.05A? — .94 + 1.33) + i(1.43A? — 3.00A + 1.82), 
we set up the matrix 

D, = (10D* — 11D* — 9D + 13) + i(14D? — 30D + 18) 
53 —84 8 36 18 56 —176 = 120 
= —9 53 —4 44} + i —30 18 36 —56]. 


Since only the first two rows will be used, the others need not be written. Starting 
with a vector whose elements are power of (—1.05+1.432)—', we get the vector se- 
quence 

(1.000, —.333 —.454i, —.095 +.3027, .169—.057%), 

(1.000000, —.332486—.453252i, —.094890 +-.3014002, .168160—.057202i), 

(1.00000000, —.33248603 —.453254133, : ), 





Thus 
1 , 
hy 9 aomrctenemnnntmnien > SARNENN s qumeiiias 
— 0.33248603 — 0.45325413% 


and A3.4=1.0522167 + 0.39596097. 
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THE TURBULENT FLOW ALONG A SEMI-INFINITE PLATE* 


‘ 


BY 
P. Y. CHOU 
National Tsing Hua University, Peiping, China 


1. Introduction. The investigation of the turbulent flow along a semi-infinite 
plate is of both theoretical and technical interest, as it leads to the magnitude of the 
skin-friction acting upon the plate at high Reynolds numbers. Interpretations of the 
velocity distribution within the tubulent boundary layer have been given by Howarth! 
according to the transport theories. But the formula that fits experiment best is 
von KaérmAn’s logarithmic law of velocity distribution which is a consequence of his 
similarity theory of turbulence. 

The purpose of the present paper is to show that von Kérman’s velocity distribu- 
tion law is consistent with the program of treating turbulent flows developed by the 
author.? For the solution of the present problem we shall adopt the following three 
conditions and simplifications. 

(1) We shall use the law of dynamical similarity down-stream along the flow. 
which has been verified by Dryden’s experiments. The mean velocity along the direc- 
tion of the undisturbed flow U which is parallel to the plane of the plate when ex- 
pressed in units of the friction velocity U, has been found by Dryden to be a func- 
tion of the dimensionless variable n = yU,/v in which y is the coordinate perpendicular 
to the plate and U, is a function of x, the distance along the plate from its leading 
edge. We now generalize this result by assuming that the mean pressure # and all the 
other velocity correlations when expressed in units of proper powers of U, are dimen- 
sionless functions of y only. 

(2) Reynolds’ equation for mean motion and the equations of double and triple 
correlations will be used. In the equations of triple correlation we shall disregard the 
quadruple correlations. The range of validity of this approximation and the possi- 
bility of its experimental verification will be discussed in Sec. 5. 

(3) We shall assume the constancy of the micro-scale of turbulence \ and of the 
squares of components of velocity fluctuation across the boundary layer. This con- 
dition is used for mathematical simplicity. It is, however, consistent with the result 
to be shown below that the shearing stress is constant across the boundary layer and 
can therefore be regarded as a condition of similarity. On the other hand, as will be 
seen presently, the variation of the mean squares of velocity fluctuation can be ac- 
counted for, though their solution will be rather complicated. This part of calculation 
is omitted due to lack of measurements. 

The logarithmic law of temperature distribution in the boundary layer of a heated 
plate can also be established by the present method of approach to the turbulence 
problem and its solution will be indicated at the end of the paper. 

2. Reynold’s equations of mean motion. Let the positive half of the xz-plane rep- 
resent the given semi-infinite plane and the leading edge of the plane coincide with 
the z-axis, the y-axis being perpendicular to the plane. Using U and V to denote the 

* Received Oct. 21, 1946. 


1S. Goldstein, Modern developments in fluid dynamics (1938) pp. 361-371. 
2 P, Y. Chou, Quart. Appl. Math. 3, 38-54 (1945). 


pipet 
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x- and y-components of the mean velocity, we find the Reynolds’ equations of mean 
motion to be 





OU aU 1 Ore 

U—+V— =— —» (2.1) 
Ox Oy p oy 

aV OV 1 @ 1 dr 

pg gw Ryley (2.2) 
Ox Oy p oy p oy 


in which both the viscous stress and the terms in components of turbulent stress 
involving partial differentiations with respect to x on the right hand sides of Eqs. 
(2.1) and (2.2) are neglected. This is permissible on account of the thinness of the 
turbulent boundary layer. 

To satisfy the equation of continuity, we introduce a stream function V as usual 


and write 


ov ov 
7= —-—--: (2.3) 


oy Ox 


| 
_ 


U = 


It is well-known that by integrating Eq. (2.1) across the boundary layer, we find 
the turbulent shearing stress on the plate. Thus,® 


; 2 0 8 
7. tae |y=0 + pu, = < f p(U, —_ U)Udy, (2.4) 
Ox 0 


where U, is the friction velocity, 6 the thickness of the layer and U, the velocity in 
the free stream outside the layer. Hence if the mean velocity U, as a function of y, 
were known, Eq. (2.4) would define the friction velocity U, as a function of x. 
Here we introduce Dryden’s experimental result that the ratio of the mean 
velocity U and U, is a function of the dimensionless variable n, where 7 is defined by 


n = yU,/v. (2.5) 
In other words, the stream function V is given by 
W = vF(n) (2.6) 
so that according to equation (2.3) 
U = U,F'(n). (2.7) 


This relation can be construed as the condition of dynamical similarity in the differ- 
ent cross-sections of the boundary layer. From Eqs. (2.3) and (2.6), we obtain 


0U, 
Vy a (2.8) 
Ox 
Let us assume that dynamical similarity also holds true for the shearing stress, 
i.e. that 
1 2 
— tay = Usa). (2.9) 
p 


3S. Goldstein, loc. cit., p. 132. 
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Substituting from Eqs. (2.7), (2.8) and (2.9) into the equation of mean motion (2.1), 
we obtain 


PS ee erg, (2.10) 


We may estimate the order of magnitude of 0U,/dx by means of von Karman’s 
coefficient of local skin friction in the following way: 


a 
8 ny (2.11) 


dx x1+2/z 

where z= KU,/U,, K =0.4. 

The derivative dz/dx thus tends to 1/x for large values of z, Hence equation (2.10) 
becomes 

e Vv 

if -_— —— = PER, (2.12) 
; K Ux 
The Reynolds number R, in Dryden’s measurements‘ is of the order of 10° and F’ 
is around 20. This shows that fj is of the order 10-* and can be neglected. In short, f, 
or the shearing stress is essentially constant across the turbulent boundary layer. 

3. Equations of double and triple correlations, the mean velocity distribution. 
We suppose that the law of dynamical similarity not only holds true for the mean 
velocity distribution, but also for the double and triple velocity correlations as well. 


Hence we write 


1 Nd | 
— —f = Usfoln), 
p 
1 F 1 rs 1 * | 
= = Tez = l ifn), oo Tyy = l rf2(n), cr tacts Tzz =e l J3(n), ( (3. 1) 
p p p 
1 2 | 
— —Ty = U,fi(n); 
p 


Wi = U;ho(n), WW, = U-hi(n), ) 
W.W8 = U;sha(n), Wy = Ushs(n), } — 


ie es | 
WW? = U-hi(n), WW? = Ushs(n). | 


Since our interest at present lies mainly in the mean velocity distribution, we 
shall only put down the equations which are pertinent for its determination. When 
written in full the component equation (i=1, k=2) for the double correlation? is 


1 aV 0U 1 OT zy 1 OT zy 
wie ats Tez i = + Tyy = U ——— — V ix 
p Ox ) p Ox p oy 


a OU = vk 
+ — wiw, + 
Ox ts) 





WW, = — O22 — — — Ys Wy. (3.3) 
y Oy d? 


4H. L. Dryden, N. A. C. A. Report No. 562 (1936). 
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In the above equation we have only retained the largest term de120U/dy among the 
sum d4.U™,n- 


Next, introduce the substitutions (2.7), (2.8), (3.1) and (3.2) into Eq. (3.3). 
Obviously, the terms involving 0V/dx and ww, in (3.3) are small. With the use of 
the relation (2.11), we find that 
1 ¥ OT zy 1 y OT zy aU, 


F'fy. (3.4) 


p Ox p oy Kxu, 
The appearance of xU, in the denominator means that this term contains the re- 
ciprocal of the Reynolds number R, as a factor when compared with the remaining 
terms, says, —T,yOU/pdy, in the equation. Hence as a good approximation of (3.3) 
we have 


- aF” + hd == lof s, (3.5) 
where /) and a are defined by 
lo = 2hv?/MU,, @= — (fe+ aeus/U,). (3.6) 


The fact that only the terms in 0U/dy and dw,w?/dy on the left-hand side of 
equation (3.3) remain can be interpreted that V is negligible and all the functions 
U, 7; and w;w,w, depend upon 7 alone and are very approximately independent of x, 
the coordinate from the leading edge of the plate. An equivalent statement is that 
the turblent boundary layer is very thin. This simplifies the problem; its mathe- 
matics becomes analogous to that for the problem of mean velocity distribution in 
a channel.® 

The component equation ({=1, k=1, 1=2) of the equations of triple correlation 
which is needed for the determination of the mean velocity distribution can be ob- 
tained by the method used in deriving equation (3.5) as: 


ou | oe ou 1 OT zy OT yy - 
3 W2rWy + wrwy = — dbei12e— — Cue + — | 2tzy — H+ Tez . (3.7) 


oy Oy dy p? oy oy 





It is obvious that the quadruple correlation term dw,w?/dy must be of the same 
order of magnitude as 7,,0Ty,/p?0y and we shall now neglect these quantities against 
the term 0U/dy(2w,w?+ e112). This is equivalent to the statement that we limit our 


solution to the region where the following inequality holds true 


| aU |. dw 9_ 


| ay (2 w.w? + bere) | > ; — r (3.8) 
7 y y 


apy2 
Wrz 








Apparently this condition can be satisfied in the middle of the layer where 0U/dy 
is large, while on its outer edge it may break down as we shall see in Sec. 5. The 
quantity Ci in (3.7) will be considered as independent of » because it varies but 
slowly with the coordinates, as has been proved previously.? The derivative 072,/dy 
is small according to the relation (2.12) and can therefore be neglected. Introducing 
expressions (2.7) and (3.2) into Eq. (3.7) and writing 


5 
C12 = 2cl 7/¥, 


5 P. Y. Chou, Quart. of Appl. Math., 3, 198-209 (1945). 
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we find, after using the inequality (3.8), 
/pve / ,3 
he = << c/F — bo1112/ 2l rT (3.9) 


Like ¢y12 we take bzi112/U? to be also a constant. 
Elimnation of the variable hz between equations (3.5) and (3.9) yields the final 
equation for the determination of means velocity distribution, 


vf ee eee (3.10) 
aF [— = lof. 5.10 
Limo 


Since the variation of f, with » across the boundary layer is small, we assume that a 
which depends upon w? defined by equation (3.6) is also a constant. There is no a 
priori experimental justification for the assumption, but it simplifies the calculation 
a great deal. As a matter of fact the determination of the behaviour of w?, w? and w? 
within the layer is not wholly impossible by means of the rest of the equations of 
double and triple correlations which have not been written out explicitly,* though the 
work involved will be rather complicated. For the same reason of simplifying calcu- 
lation we take J») to be constant, namely, A to be independent of 7 and leave out the 
equation of vorticity decay.? 
Integrating the equation (3.10) twice, we find 


a 
n = —F’ + A(eoeF'le — 1), (3.11) 


oj 4 


where A is a constant of integration. The other arbitrary constant has been deter- 
mined by the condition on the surface of the plate: 7=0, F’=0. Another constant 
of the remaining three in equation (3.11) can be disposed of, as a consequence of the 
linear dependence of the mean velocity U upon 7 near the plate, by requiring that it 
become »=F’, when 7 is small. On the other hand for large values of 7 the mean 
velocity U varies with y logarithmically, confirming von Karman’s result. For this 
latter reason we take the factor of F’ in the exponential function to be 0.4, when we 
compare formula (3.11) with Dryden’s measurements.‘ The last constant A is de- 
termined by passing the theoretical curve through the mean experimental value 
of F’=18.2 at logio 7 =2.2. The result can be summarized in the equation 


L 
n= rr + 0.1024(e9-47/Ur — 1), (3.12) 

The velocity distribution given above resembles that within the channel in its 
functional dependence upon 7, though the values of the corresponding constants ap- 
pearing in the two cases are different as is to be expected. Comparison between 
theory and experiment is given in Fig. 1. Their agreement is well-known to be 
satisfactory. 

* They will be the same as Eqs. (2.3), (2.5), (2.6), (3.1), (3.3) and (3.4) in the paper quoted in Foot- 
note 5. 
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Fic. 1 


Velocity distribution in a turbulent boundary layer. 





Theory: - ——— 
Dryden’s experimental results: 

R, = 335,000, (] R.=507,000, 

< R.,=788,000, A R =983,000, 


x being measured from the transition point. 


4. The coefficient of local skin-friction. The explicit relation between the coeffi- 
cient of local skin-friction cy, which is defined by 
2 : 
cy = 70/2pU1 = U,/3U,, (4.1) 
and the Reynolds number R,= U,x/v can now be determined by substituting the 
mean velocity distribution (3.12) into Eq. (2.4) and integrating. The thickness of the 
boundary layer y=64 is found by setting U=U, in Eq. (3.12). Retaining the large 
quantities in the final result after integration, we find 
dR, 
~ = 1.6e*(z? — 2s + 2), (4.2) 
dz 
where 2 has been defined in Eq. (2.11) as 0.4U,/U,. Integrating (4.2) again gives 
R, = 1.6e*(z? — 42 + 6) — 6. (4.3) 
The constant of integration is determined by x =0, when z=0, under the assumption 
that the flow is turbulent from the leading edge. 
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For large z we keep only the 2? term on the right member of equation (4.3) and 


a 


the result can be put in the form 
cp” = 1.19 + 4.08 logis (RC). (4.4) 


The two constants appearing above according to Kempf’s experiments are 1.7 and 
4.15 respectively.! The form of R, given in (4.3) also justifies our use of the relation 
(2.11) which leads to the very approximate constancy of the shearing stress across the 
boundary layer as has been demonstrated by relation (2.12). 

5. Some further experimental verifications and validity of the solution. The 
present analysis shows that, according to (2.12), the shearing stress is constant across 
the turbulent boundary layer, if the Reynolds number of the flow is sufficiently high. 
This result can be subject to direct experimental proof, since 7,, at different points 
of the fluid is now measurable by the hot-wire technique. On the other hand to 
establish the validity of the inequality (3.8) experimentally is probably not simple. 
There is a greater chance of its fulfillment, if w? and w? are also constant across the 
layer. This constancy of w2 and w? can be interpreted as the consequence of the 
condition of “similarity” mentioned before. On the other hand, even if this require- 
ment of similarity were satisfied, it has its limitations and we can determine the 
validity of the present solution by considering the values of the constants appearing 
in equations (3.11) and (3.12). 

Since f, is independent of 7, according to its definition (2.9) and relation (2.4) it 
must be —1. From Eggs. (3.11) and (3.12) we find 


a/lg = 0.959, lo/c = — 0.4. +3 1) 





The quantity /» is determinable by using Eq. (3.6) with the micro-scale of turbulence 
\, which comes from correlation measurements and k, a constant, which is present 
in all the three equations for w?, w?, and w? and is also observable. Thus both a and c 
are known after /) is fixed. 

The factor 1/F’’ in equation (3.9) for hz varies according to (3.12) between unity, 
when U/U,=F’=0 on the plate and about 120 when F’ =20. Hence the part of he 
variable with 7 can take the extreme values of —c and —120c. According to its 
definition in (3.2), the value of hk, can not be very large. This means that there must 
be an upper limit of F’, namely, the outer limit of the boundary layer within which 
the relation (3.9) is valid. Again if triple velocity correlations were measurable by 
the hot wire technique, to test this relation directly within its region of validity should 
be an interesting problem. 

The foregoing analysis shows the important réle played by the triple velocity cor- 
relation within the turbulent boundary layer. This can be interpreted as indicating 
that turbulence is generated in this layer and transported away from it. The triple 
velocity correlation is a measure of this transport; hence the set of equations govern- 
ing its behaviour is needed for the interpretation of the mean velocity distribution. 
Indeed, if this interpretation were correct, it is easy to understand that only the 
equations of mean motion and of double correlation with the omission of the triple 
correlation terms in the latter are insufficient to explain the velocity distribution 
within the boundary layer, for this simplified set of the equations of double correla- 
tion only tells us the balance between the production of turbulent energy by de- 
formation of the mean flow and the decay of turbulence. Although this picture 
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seems to give a good account of the velocity distributions in the center of the channel 
and in free turbulence like jets and wakes, it is definitely insufficient for the phe- 
nomena dealt with in the present paper. 

The present paper is obviously incomplete in the sense that theoretical computa- 
tions of the mean squares of turbulent fluctuations are not given. It does show, how- 
ever, that a certain type of function for the mean velocity distribution is necessary. 
Though refinements of calculation and measurements may change this result to some 
extent, the main feature of the function will probably remain the same. This point 
has been borne out by the corresponding calculation of the velocity distribution in 
the channel under the assumption of different turbulent levels.® 

The temperature distribution within the boundary layer of a heated plate can 
be calculated by means of the equations of mean temperature and of the double and 
triple velocity and temperature correlations. We obtain an equation for the mean 
temperature distribution analogous to (3.11), and it agrees with the measurements 
of Elids.° Obviously a limitation on the validity of the solution similar to (3.8) is 
also present. 


6 F. Elias, Zeits. angew. Math. Mech., 10, 1-14 (1930). 











—NOTES— 
A CONTRIBUTION TO THE METHOD OF LEAST SQUARES* 


By M. HERZBERGER aAnp R. H. MORRIS (Eastman Kodak Company) 


In trying to calculate the image figures of an optical problem the authors were led 
to consider a problem of least squares which proved unusually difficult, since the de- 
terminant of the normal equations to be solved had a very small value. The following 
attempt to make the problem a geometrical one should be of value in other applica- 
tions, and is presented here for what it is worth. 

Most problems of least squares can be reduced to the following: Suppose that we 
are given a sequence of k known functions fi, - - - , f; of one or more variables. Our 
problem is to determine k constants, C;, so that an unknown function, F, whose values 
are known for n > values of the coordinates will be approximated as closely as possi- 
ble by a linear combination of the functions f;: 


F=Cifit--++Cife (1) 


Since F is given for only nm >k values, we can consider these values as the coordinates 
of a vector 6 in an n-dimensional space. The values of the f; for the ” values of the 
coordinates may also be considered as the coordinates of k vectors d; in the n-dimen- 
sional space. It is a simple geometrical problem to find the best approximation of the 
vector 5 by a linear combination of the vectors d;. We merely project 5 in to the space 


of the vectors 4;, i.e., we minimize 
k 2 
(3  : Ca). (2) 


i=l 
This leads to the k normal equations for the constant C; 


bis C;a;d) = bay, (3) 


which must be solved. 

We can simplify the problem considerably if we order the vectors in a suitable 
way, which in this case will be the transformation of (3) into a form which contains 
only unit vectors, i.e., the projection of our problem on the unit sphere. Let 





G;=a5;, 6 = pi, (4) 
where sj =/2=1. Equation (3) becomes 
Dd visisn = ts. (S) 
A solution of (5) which yields y; will also yield C; for 
C= =, "= ne (6) 
ai B 


From (4) and (5) we see that the diagonal terms of the matrix of (5) are all 1, whereas 
the members not on the diagonal are smaller than 1. 





* Received Nov. 1, 1946. Communication No. 1123 from the Kodak Research Laboratories. 
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The senior author now suggests an ordering of the vectors in such a way that the 
principal subdeterminants* have a maximum value (which will always be less than 1, 
since k independent edges of the parallelopiped formed by the unit vectors §; lie inside 
the unit sphere). If the determinant of the problem is small, this means that the vol- 
ume of the parallelopiped is small or that the vectors 4; are not “too independent,” 
i.e., that one or more of the vectors can be approximated by a linear combination of 
the others. 

If we denote by a;; the quantity, s,s;, by c& the quantity ¢s;, we may write (5) in 


the form 
(0) 
aiY1 + a12Y2 foot + QurYik = C1 , 


(7) 
(0) 
GurvYi + Aneye2 + +++ + Oeeve = Cx 
The method used will be this: We choose an arbitrary number, say r, where r<k, 
and find the r-rowed principal subdeterminant which has the largest value. We set 
the (k—r) variables not involved equal to zero, and solve the r equations for the r 
unknown variables, obtaining the values y;1=Yi, * + * Y¥r1,¥31=0, j>r. Substitut- 


ing these for the variables y1, - - + ¥;, in the k—r remaining equations, we have 
(1) 
Q rt1Yr41 +ee tay =HCi , 
a) 
Qk r+ 1YVr4+1 + Sn + QE KYk = Ck, (8) 


where cf? =c —}°4_,a:,77;,. At this point we have two possibilities. Either the 
constants c}” are very near to zero, in which case we have that the solution of (8) is 
given by y;=0,7>r7, so that a reasonably accurate solution of (7) and thus of (5) 
is given by 

Vii = Via (¢=1,2,---,7); Y¥1=0 (@=r+1,---, R). (9) 


In the event that they are not sufficiently small, we take the largest (r+1)-rowed 
principal subdeterminant of (7), equate the remaining (k—r—1) variables not in- 
volved to zero, and proceed as before, obtaining (k —r—1) equations in (k —r—1) un- 


knowns: 
k—r—-1 (2) 
DS vara = 6 ’ (¢=1,2,---,&-—9r=— 1), (10) 
A=1 


where the meaning of the terms 4 and ¢” is clear. Should the constants c) be suffi- 
ciently small, the approximate solution is given by 


Vi = Yi (§ = 1,2,---,¢7+ 1); vi = 0, (¢=r+2,---, k). (11) 


One must keep in mind that as larger subdeterminants are taken, their values get 
smaller, and that, of course, as r gets nearer k the c™ become smaller, and vanish 
when r=k. There will usually be a place somewhere in this process, however, where 


sufficiently good accuracy will be achieved. If greater accuracy is desired, one must 


* Subdeterminants whose principal diagonals form part of the principal diagonal of the determinant 


of the system (5). 
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attempt one of the ordinary approximation methods, in particular, an iteration pro- 
cedure. 

We present as an example a set of six equations with six unknowns. Since the cal- 
culation of the maximum subdeterminant is very laborious, we have replaced it by 
the less laborious, but theoretically less satisfying, method of ordering the vector so 
that the sum of the squares of the elements outside the principal subdeterminant is 
as small as possible. 

The system is represented by the following matrix. 

I 


6.35083641  3.38272218 1.83936868  1.02023449 0.57681127 0.33173876 2 .27198287 
1.83271729 1.01325021 0.57097916 0.32757965 0.19093008 1 .26347780 

0.56915637 0.32551164 0.18928551 0.11166820 0.71717786 

0.18871015 0.11107961 0.06624255 0.41473305 

0.06609478  0.03979216 0.24391707 

0.02415650 0.14552476 


This matrix is reduced to unit diagonal and reordered as described above (in this ex- 
ample, the ordering remains the same). The elements of the resulting matrix are the 
light-face elements of the following matrix. 





Il 
| 1.00000000 0.99152377 0.96747004 0.93193774 0.89029657 0.84696167 0.90154952 0.90154952 
| 490000000  —- 0.99209478 0.97090085 0.94120796 0.90742269 0.93329702 
—0.99152377 | 0.01688061 0.03282524 0.04686243 0.05845775 0.06764006 0.03938924  2.33340146 
ii | 1,00000000 0.99323746 0.97592846 0.95235029 0.95062903 
—0.96747004 —1.94455295 | 0.00017130 0.00048914 0.00091901 0.00141057 0.00181242 10.58038529 
ae | 1.00000000 0.99461070 0.98112063 0.95470899 
—0.93193774 —2.77610999 —2.85545826 | 0.00000007 ' 
| 1.00000000 0.99585694 0.94876491 
—0.89029657  —3.46301170 —5.36491535 
- 1.00000000 0.93631051 
—0.84696167 —4.00696776 —8.23450088 | 


The bold-face elements are found as follows: The first bold-face row contains the 
determinants of the second order found from the first and second rows. Any other 
bold-face elements above the stairway dividing line are formed by adding to the ele- 
ment immediately above it the product of the rest of the column above it and the row 
to its left, e.g., 0.99323748 —0.93193775 X0.96747003 — 0.04686244 X 1.94455183 
= ().00048919. The first column to the left of the stairway dividing is the negative of 
first row, the second column is the negative of the first bold-face row divided by its 
first number, and so on. The process is continued until one of the diagonal elements 
becomes very small. The solution of the equations, omitting the last and successive 
variables, is found in the usual manner from the bold-face column at the extreme 
right. This column is the ratio of the last number of each row to the first number of 
each row. 

We find in our case’ 


v1 = 8.75144919, -y2 = — 18.24071797, ys = 10.58038529. 
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The three remaining unknowns are found by distributing the residuals among all six 

equations by least squares. When this is done, we find for our original unknowns: 
C, = 3.4727, Cz = — 13.4739, C; = 14.0244, 

— 0.0003, Cs = 0.0002, C. = 0.0001. 


C4 


This method has been applied successfully to least-squares solutions in geometrical 
optics and to the colorimetric problem of finding polynomial reflection curves which 
will yield a prescribed set of tristimulus values under a given illuminant and fit other 
prescribed conditions. 

In the latter problem, the small unknowns must be neglected prior to symmetriza- 
tion and reduction to unit diagonal. 


NEW FORMULATIONS OF THE EQUATIONS FOR 
COMPRESSIBLE FLOW* 


3y B. L. HICKS (Ballistic Research Laboratories, Aberdeen Proving Ground) 
P. E. GUENTHER (Case School of Applied Science) anv R. H. WASSERMAN (University of Chicago) 


Introduction. A prominent aerodynamic effect of combustion in a moving gas 
stream is an alteration of the flow pattern owing to heat release within the fluid. 
This alteration occurs not only in the immediate neighborhood of heat sources but 
also downstream where the entropy and stagnation temperature vary from one 
streamline to another. As a background for combustion research, appropriate de- 
scriptions of these altered flow patterns have been investigated. This paper considers 
the downstream patterns, which are restricted to be the adiabatic and steady flows 
of an inviscid fluid. In a second paper,’ diabatic (i.e., non-adiabatic) flows will be 
discussed. 

Since one-dimensional flow theory?:* can be considerably condensed by use of the 
local Mach number M, it was natural to seek a corresponding condensation with the 
help of the Mach vector 


M = MV/V 
and the Crocco vector 
WwW = V/V: 


in which V, is the variable limiting velocity at each point of the fluid. The introduc- 
tion of Mach and Crocco vectors into the compressible flow equations sufficiently 
simplified or altered their form that a number of further investigations were sug- 
gested including those of diabatic flow. 


* Received Feb. 12, 1947. This paper is a revised report of theoretical work performed by the authors 
at the Cleveland Laboratory of the National Advisory Committee for Aeronautics in 1943-45. 

1B. L. Hicks, Diabatic flow of a compressible fluid, submitted to Quarterly of Applied Mathematics. 

? Neil P. Bailey, The thermodynamics of air at high velocities, Journ. Aero. Sci. 11, 227-238 (1944). 

3B. L. Hicks, D. J. Montgomery, and R. H. Wasserman, The one dimensional theory of steady com- 
pressible fluid flow in ducts with friction and heat addition, NACA TN, 1947. 
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1. Development of the new equations. The basis for our development is the fol- 
lowing set of equations: 


cpV-VT + pv-V = 0, (1.1) 
pV-VV+ Vp = 0,7 (1.2) 
V-pV = 0, (1.3) 

p = RT. (1.4) 


It is appropriate to include (1.1), an expression of energy conservation, because the 
entropy S and stagnation temperature 7; are not constant throughout the field of 
flow. The Bernoulli relation, 

Cpls = Cy + $V’, (1.5) 


then defines the stagnation temperature. In terms of M and W, equation (1.5) be- 
comes 
(T7,/T) = 1+ (vy — 1)M?/2 = (1 — W?). (1.6) 
The Mach vector equations are derived from (1.1), (1.2), (1.3), (1.4) by the sub- 
stitution V=*/yRT M and elimination of p and T. The new equation of motion re- 
sulting is 


y¥-—1 1 
M-VM — —— MV-M + — V log p = 0, (1.7) 
eT Y 
and the new equation of continuity 
Y— 1 —(y+1)/2(y—-1) 
v-(1 +— u) M = 0. (1.8) 


The equations of motion and continuity in terms of the Crocco vector are derived 
similarly, the initial transformations being V= V,W and T= V7(1—W?)/2c,. It is in 
addition necessary to use the fact that JT; and therefore also V; are constant on stream- 
lines for adiabatic, inviscid flow (cf., for example, Vaszonyi’s proof*). The Crocco 
vector equations are then 
i aig ” 
W-vw + — (1 — W2)V log p = 0, (1.9) 


V-(1— W2)" Ww = 0. (1.10) 


2. General comments. The local behavior of a compressible fluid is characterized 
by its Mach number. Therefore the vectors M and W, from whose magnitudes M can 
be determined, supply more useful descriptions of a flow than the vector V from whose 
magnitude M cannot be computed without knowledge of V;. The equations possess 
equal generality in all three languages but are simpler in M and W language for these 
representations have permitted elimination of all but one vector variable from the 
equation of continuity and elimination of all but this vector and p from the equation 
of motion. Such reduction in the number of independent variables has heretofore been 
accomplished only for iso-energetic flow. 


4 A. Vazsonyi, On rotational gas flows, Quarterly Appl. Math. 3, 29-37 (1945). 
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Rotational flows are most aptly described by the W language; the added term 
MV -M in the M equation of motion brings with it no advantage. (The W language is 
still the best when diabatic rotational flows are considered.) However when either M 
or W is irrotational, the physical characteristics of the resulting flows make it ap- 
propriate to retain each language. 

3. Geometry of stream tubes and transition through the speed of sound. In the 
one-dimensional approximation, isentropic compressible flows have minimum flow 
area at sonic velocity. The same property has been suggested by Prandtl® for irrota- 
tional three-dimensional flows. We shall prove that this property is valid for all con- 
tinuous, steady adiabatic flows when an appropriate description of stream-tube area 


is introduced. 





Fic. 1. 


Let B be any continuous and solenoidal vector parallel to M. Consider any two 
plane sections AA;, AA» of the stream-tube 2 which are normal to the streamline L 
contained in 2. Since Jaa,B-dA, = fx4,B-dAs=fx4B,dA a position (m) can be found 
for each section such that [B,(m)AA hh = [B,(m)AA |]. or 

i B,(m) + — = 0 (3.1) 
— log B,(m —_ — = ‘ 
Os AA as 
in which 0/0s denotes differentiation along L. As AA-0, B,(m)—>B, the magnitude of 
B on L. Also if s=B/B is the unit vector in the direction of flow 


1 odAA 


im — 3.2) 
44-0 AA Os ( 


V's = 





The fractional rate of change of stream-tube area with respect to streamline arc is 

therefore measured by V-s for flow in both two and three dimensions. Similar relations 

derived geometrically have been used by von Ka4rman® for two-dimensional flows. 
The continuity equation in M language can be rewritten in terms of V-s as 


y-1_ \" aM 
(1 ane ae u) lll kt hi (3.3) 
; s 


For V-s>0, |1—M| cannot decrease but for V-s <0, |1—M| must decrease in the di- 


® L. Prandtl, General considerations on the flow of compressible fluids, NACA TM No. 805, 1936. 
6 Th. von Karman, Compressibility effects in aerodynamics, Jour. Aero. Sci. 8, 337-356 (1941). 
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rection of flow i.e. sonic velocity is approached. As M-—>1 however, 0M/ds diverges 
like —(y+1)(V-s)/4(1—M) and 0/0s diverges similarly. Consequently if M is to 
pass through the value unity continuously, V-s must be zero there. Consideration 
of the signs shows that the indicated extremum corresponds to a minimum of stream- 
tube area at M=1. 

4. Thermodynamic considerations. The equation 


VX (V XV) =¢,VT; — TVS (4.1) 


which is similar to Vazsonyi’s,* contains the thermodynamic implications of the V 
equations. (Less general equations than (4.1) have also been derived by Crocco,’ 
Tollmien* and Emmons.*) The corresponding equation in terms of W and the stagna- 
tion pressure ; is easily obtained from (1.9) 


wx (VxXW) = <— (1 — WV log p, (4.2) 
showing that there is one fundamental thermodynamic quantity in adiabatic flow, 
namely ;, and not two separate ones S and 7;. Thus for a given p,=/ distribution, 
the same W (or M but not V) field obtains no matter what the variation of Sand T, 
so long as they combine in the form [T,exp(—S/c,) |" to give p/. There seems to 
be no corresponding theorem for equation (4.1). 

Differentiation of (4.1) yields a differential form of Bjerknes’ theorem® 


vx [Vx(VxV)] = —-—VITx VS. (4.3) 
Differentiation of (4.2) results in Crocco’s equation’ containing W only: 
[ee 0 


bod 4.4 
1— W? — 


Equations (1.10), (4.3) and (4.4) were originally derived for isoenergetic flow, a 
restriction now seen to be unnecessary. 

5. Irrotational flows. We now consider in turn the consequences of irrotationality 
of M, W and V fields. If M=Voy, an integrability condition 


Vom = f(¢m) (5.1) 


must be satisfied if the M equation of motion (1.7) is to be integrable. (A similar 
equation for the stream function occurs in theory of rotational incompressible flow.) ". 
Since the potential function must also satisfy the continuity equation (1.8) it is to be 
expected that the manifold of permissible functions gy will be somewhat restricted. 
This is illustrated by the fact that p exp (yM?/2) isa function of gy alone, or that 
M, om and p; are related by the expression 

7 L. Crocco, Eine neue Stromfunktion fir die Erforschung der Bewegung der Gase mit Rotation, ZaMM 
17, 1-7 (1937). 

®W. Tollmien, Ein Wirbelsats fur stationdre isoenergetische Gasstrimungen, Luftfahrtforschung 19, 
145-147 (1942), British R.T.P. Trans. No. 1744, Ministry Aircraft Prod. 

9° H.W. Emmons, The numerical solution of compressible fluid flow problems, NACA TN No. 932, 1944. 

10M, Lagally, Ideale Flissigkeiten, Handbuch der Physik, Julius Springer (Berlin), vol. 7, ch. I, 


Art. 19, p. 29, and Art. 32, p. 49. 
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1. 1 

=(M) = &(¢m) ae ee log Pu (5.2) 
7 


in which 


— 1 a 1 —1 
"(M) = —— — u(1 + — u) dM, 


a 


—1 
P(yu) = 1 J Kewdew. (5.3) 


(In uniplanar flow, explicit formulae for potential and stream functions can be derived 
for both irrotational M and some irrotational W flows. We have not ascertained 
whether these formulae represent more than the well-known elementary radial and 
vortex flows. ) 

If W=Vow, then Vp;=0, a case discussed by Vazsonyi‘ in the less appropriate V 
language. Flows with constant stagnation pressure but variable entropy and stagna- 
tion temperature from one streamline to another may be of immediate practical in- 
terest. For example the effect of variable chamber gas temperature upon thrust co- 
efficient of a rocket with “perfect” nozzle can be computed. This calculation is possi- 
ble because irrotational V flow is the special case of irrotational W flow which occurs 
when not only ~,; but also S and therefore T; are constant throughout the flow 
(VS=VT,=0 is implied by V=Vgy). Accordingly the same partial differential equa- 
tion (derived from the continuity equation (1.10)) is satisfied by both gw and its 


isentropic form ¢gy, 


: 2 IPw 271 3? w + : OGw 0 w 3 'w 
3 ee = | ce 


inh y —1\0x; Ox? ¥ —1 jer Ox; Ox; Ox,0x; 


LINEARIZATION OF SOLUTIONS IN SUPERSONIC FLOW* 
3y JOHN W. TUKEY (Bell Telephone Laboratories and Princeton University) 


1. Introduction. The equations governing flow at supersonic speeds are believed 
to be well known, but the difficulties of calculating exact solutions are so great that 
approximate solutions are the aim of the present and foreseeable future. Two main 
approaches to such approximate solutions are commonly considered: 

(a) the calculation of numerical approximations by high-speed calculators, 

(b) simplification of the equations and explicit solution of the simplified equations. 
Of course the first route requires the high-speed calculating machines which are now 
in sight, but not yet available. 

It is the purpose of this note to propose and exemplify a third approach which 
may partly replace (b) and frequently supplement (a), namely 
(c) simplification of the dependence of the answers on one or more parameters. 

2. Is linearization of solutions mathematically justified? The first objection which 
many would raise against “linearizing the solutions” as a substitute for “linearizing 


* Received Feb. 15, 1947. 








362 NOTES [Vol. V, No. 3 
the equations” is its mathematical justification. The justification of the former would 
seem to be at least as good as the latter. In order to make this clear, let us look at the 
situation in more detail. 

A problem is believed to be adequately stated by a complex set of equations, which 
we may symbolize by 


F(x, ) = 0, (1) 


where \ is a parameter. The classical procedure is to replace the complex set of equa- 
tions by a simpler set, frequently composed of linear equations, which we may write 


L(x, ) = 0, (2) 
and then to solve these equations exactly to obtain 
x = g(A). 


We must interpret this solution as an approximate solution of the original equations. 
What justification might we have? The only possibilities seem to be 
(a) a theorem that the solutions of (2) are close to the solutions of (1), 
(b) empirical evidence that, in particular cases, the solutions of (2) are close to the 
solutions of (1), 
(c) mathematical intuition, 
(d) physical intuition. 
The author knows of no case in supersonic flow where tke very desirable justification 
of type (a) exists. 
On the other hand, suppose that (1) can be solved for, say, two values of \, so 
that we know x; and x2, where 


F(x1, \1) = 0 = F(xz2, do), (3) 

and suppose further that we believe that we can approximate the dependence of x 
on \ over a certain range by one of, say, a two-parameter family of functions, 

x(A) = h(d; A, B). (4) 


We can then use the two known cases to determine the parameters and find an ap- 
proximate solution 


wal 


x = h(d; Ao, Bo) ‘4 


of (1). What justification might we have? The only possibilities seem to be 

(a) a theorem that the solution of (1) depends on d in a way close to one of the func- 
tions (4), 

(b) empirical evidence that, in particular cases, the solutions of (1) do behave like 
one of the functions (4), 

(c) mathematical intuition, 

(d) physical intuition. 

The author knows of no case in supersonic flow where the very desirable justification 

of type (a) exists. 
Thus the two ways of “linearizing” seem equally justified. The first appeals to 

some because the solution used is exactly a solution of (2). There is no reason, how- 
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ever, why this should be more esthetic than the fact that the other approximation 
is in exactly the form 4. 

3. A linearization for the region of attacked shock waves. In the range of super- 
sonic speeds safely away from sonic or hypersonic speeds, there seems to be some evi- 
dence that a linearization of the form 


where g, the dynamic pressure, is given by 
7 
q = 2 M*p, (6) 


applies approximately to many averaged quantities, such as lifts, drags, etc., for well 
streamlined bodies. Some empirical evidence is presented in the next section. It is 
believed that, when the two points needed to determine A and B are known from 
theory or experiment, the use of this linearization will frequently give results prefer- 
able to those obtained by solving a linearized problem. 

4. The flow around an axially-symmetrical cone. The theoretical problem of flow 
around a cone was solved by Taylor and Maccoll! who computed the drag or pressure 
coefficient for several cases. These calculations were refined and extended by Maccoll.” 

A linearized theory was previously presented by v. Karman and Moore.*® 

The present paragraph is concerned with the empirical linearization 


0.096\ / 6, \1-69 
C, = (0.083 + —— )(<5) : (@, = half angle of cone). (*) 
M? 10° 

The form and constants of this approximation are purely empirical, and are princi- 
pally based on Maccoll’s theoretical values. 

Figure 4 on p. 309 of the paper by Taylor and Maccoll? shows that the linearized 
equations give rise to errors of 10 to 25 per cent over substantial portions of the 
range, while the following table shows that the empirical formula gives agreement 
to better than 5 per cent over a range of half angles from 10° to 55° (Maccoll estimates 
57.4° as the maximum half angle for conical flow) and over a range of Mach numbers 
for each angle extending down to within about 0.1 of the minimum Mach number for 
conical flow (attached shock-wave). 

The fit is not perfect, and it is clear that the empirical formula cannot replace the 
exact solution for all purposes. However, as a means of interpolating between known 
values, it clearly serves a useful purpose. 

Summary. There is theoretical and empirical justification for “linearizing the 
answer” rather than “linearizing the equations” in some fields of supersonic flow. 


1G. I. Taylor and J. W. Maccoll, The air pressure on a cone moving at high speeds, Proc. Roy. Soc. 
London (A) 139, 278-311 (1933). 

2 J. W. Maccoll, The conical shock wave formed by a cone moving at a high speed, Proc. Roy. Soc. London 
(A) 159, 459-472 (1937). 

3 Th. von Kérm4n and N. B. Moore, Resistance of slender bodies movit'g with supersonic velocities, with 
special reference to projectiles, Tr. Amer. Soc. Mech. Eng. 54, 303-310 (1932); for further discussion see 
Taylor and Maccoll, loc. cit., p. 30. 
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TABLE 1. Supersonic Pressure Coefficient of Unyawed cone 








Semiangle of Cone Mach Pressure Coefficient! Per Cent 
(Min. Mach No.?*) Number Quadrature*® Empirical* Deviation 

10° 1.091 1565 .161 +5 
(1.04) 1.403 .130 131 +1 
1.81 .108 442 4 
2.39 .096 .100 +4 
3.33 .090 .092 +2 
5.46 .076 .086 +13 
20° 1.314 .446 -447 +0 
(1.18) 1.555 .374 394 +5 
1.65 .356 382 +7 
243 .318 333 +5 
2.86 .296 306 3 
4.17 268 286 +7 
9.74 -256 271 +6 
30° 1.514 .814 .800 —2 
(1.46) 1.639 .730 .760 +4 
2.05 .636 .678 +7 
pe | .598 614 +3 
3.16 .570 .593 +4 
3.85 .556 572 +3 
7.84 .524 542 +3 
40° 2.03 1.098 1.106 +1 
(1.95) yy 1.036 1.067 +3 
2.84 0.948 0.988 +4 
4.20 0.904 0.921 +0 

5.86 0.890 0.893 
45° 2.38 1.316 1.270 —4 
(2.3) a i 1.182 1.219 +4 
3.76 1.104 1.140 +4 
4.89 1.084 1.105 +2 
8.01 1.068 1.073 +1 
50° 3.45 1.440 1.407 —2 
(3.1) 3.33 1.374 1.391 +1 
3.81 1.330 1.360 +2 
4.75 1.296 1.325 +2 
7.01 1.266 1.289 +2 
55° 5.59 1.552 1.535 —1 
(5.5) 6.38 1.516 1.522 0 
7.56 1.496 1.510 +1 








1 Given in American units—Maccoll’s values have been doubled. 

2 Minimum for conical flow and attached shock wave. 

3 Taken from Maccoll, Proc. Roy. Soc. 159 (1937) pp. 461 and 469, also 139 (1933) p. 285. 
* Calculated by formula (*) above. 

5 This value is undoubtedly too low by 10 to 20 per cent. 
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The linearization 
(drag, lift, etc.) = Ap + Bq 


is proposed for use in the region of attached shock waves. 


BOOK REVIEWS 


Tables of spherical Bessel functions. Prepared by the Mathematical Tables Project, 
National Bureau of Standards. Vol. 1. Columbia University Press, New York, 


1947. xxviii+375 pp. $7.50. 


In this volume the functions \//2xJ,(x) are tabulated for the half-odd-integer values of », 
—27/2Sv527/2. For 0<x310.0, the functions are tabulated to eight or more significant figures at in- 
tervals of .010. For 10<x<25, the tables contain seven significant figures and the argument interval 
is .10. A table of the function x? tl24 /¢ /2xJ,(x) is also included for argument intervals of .01. The range 
of argument is: for y= —4, +5/2, +7/2:0sx381; for »= +3/2:0<5x8.50; for y= +9/2, +11/2,---, 
+15/2:0<x<1.5; for v= +17/2, -- +, 23/2:05x52; and for y= +25/2, +27/2:05x%82.5. 

G. F. CARRIER 


Introduction to aerodynamics of a compressible fluid. By Hans Wolfgang Liepmann and 
Allen E. Puckett. Galcit Aeronautical Series. John Wiley & Sons, Inc., New York, 
1947. ix+262 pp. $4.00. 

The better known works on the aerodynamics of compressible fluids are 1) Ackeret’s and 2) Buse- 
mann's articles in the Handbuch der Physik and the Handbuch der Experimentalphysik; 3) Bateman’s 
article in the Report of the Committee on Hydrodynamics of the National Research Council; 4) Taylor 
and Maccoll’s article in Durand’s Aerodynamic Theory; 5) papers collected in the Proceedings of the 
Fifth Volta Congress and 6) Sauer’s rather recent book. The first three, although excellent summaries 
of information up to the dates of their publication are now outdated. The fourth is masterly in its choice 
of subjects and exposition but may be too difficult for a beginner in the subject due to its briefness. The 
fifth contains truly a wealth of good papers covering the whole field, but unfortunately it is out of print. 
The sixth and last one is written by a mathematician from the point of view of a mathematician. 

Therefore for aeronautical engineers now entering upon a study of aerodynamics of compressible 
fluids, there was a great need for an introductory book with a clear exposition of the physical aspects and 
the recent advancement of the subject. This need is believed to be admirably satisfied by Liepmann and 
Puckett’s book. Of course, as the title indicates, the book is not a compendium of all available data on 
compressible flows, and for engineers who are desperately looking for design information, this book will 
be a disappointment. However, for those who are interested in a clear understanding of the phenomena of 
compressible flow, the outstanding feature of this book is the constant emphasis on the interplay between 
theory and physical facts. The authors always try to present a description before going into a calculation. 
Thus the gist of a problem is not lost in a maze of equations. It also is evident that the authors are always 
aware of the limitations imposed by simplifying assumptions which are, of course, necessary in order to 
make the problem tractable, and such limitations are always pointed out in due time. The reader is thus 
saved from the various pitfalls in applying the results of a theory. The authors’ research efforts in com- 
pressible flow are also evident in the book because not only are known facts expounded, but difficult 
questions in such subjects of current interest as condensation shock, Mach reflection, influence of viscosity 
in transonic flows are also discussed. 

However there are occasions in the book where the authors leave an inquisitive mind unsatisfied. For 
instance, the treatment in Part II on two- and three-dimensional motions is uneven in its pace.—While 
great details of ccmputations are given in deriving the basic perturbation equations, no method of solving 
these equations or example is given, so that reader could have a more concrete understanding of the prin- 
ciples involved. Then for condensation shocks, there is the question of the relative importance of dust 
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nucleii and nucleii formed by statistical fluctuations from thermodynamical equilibrium, and the “theo- 
retical” possibility of the unobserved condensation shock to subsonic velocities. The authors’ discussion 
on these questions are vague and indecisive. The form of expanison of potential into a series of thickness 
ratio of the body in Chapter 8 is, of course, incorrect for any three-dimensional axially symmetric flow 
and for any two-dimensional flow with a stagnation point. However the reviewer does not consider these 
statements as criticisms of this excellent book.—If the book were without such defects, then it would be 
a book by a master. But such definitive treatise on any one subject can be expected to appear, perhaps, 
only once in fifty years. 


H. S. TstEn 


The theory of mathematical machines. By Francis J. Murray. King’s Crown Press, 

New York, 1947. vii+116 pp. $3.00. 

This book discusses the operation of many basic units used in mathematical machines. The interest 
is in the basic devices themselves rather than in the problems arising when they are interconnected. 

The four chapters of the book are entitled Digital Machines, Continuous Operators, The Solution of 
Equations, and Mathematical Instruments. Each of these is a series of discussions of various devices. 
A typical discussion contains a description of the physical appearance together with a brief analysis of the 
operation. Both mechanical and electrical schemes are discussed under the assumption that the reader has 
a reasonable acquaintance with the operation of simple machine components and with elementary electri- 
cal, and particularly electronic, principles. 

The subject matter is a collection of material which for the most part is widely scattered in the litera- 
ture of various fields. For this reason, the book would be extremely interesting to anyone wishing to ob- 
tain an acquaintance with the subject. The treatment of analogue devices using electrical networks and 
mechanical devices is fairly complete and quite interesting—particularly that part dealing with the au- 
thor’s own work. On the other hand, more space might have been devoted to digital techniques and less 
to the discussion of engineering design features which are occasionally discussed at some considerable 
length. 
The book is issued by the King’s Crown Press in loose leaf form “—as submitted by the author, with- 
out the usual editorial attention of the Columbia University Press,” in an effort to minimize the cost to 
the buyer. It is this reviewer's opinion that further editorial work would have been well warranted since 
there are many minor errors and vague statements, and the lack of numbered references to diagrams 


makes it difficult at times to follow the technical arguments. 
J. A. KRUMHANSL 


Methods of mathematical physics. By Harold Jeffries and Bertha S. Jeffries. Cam- 
bridge: At the University Press, 1946. vii+679 pp. $15.00. 

“Methods of Mathematical Physics” is a very comprehensive work dealing with the mathematical 
techniques of classical and modern physics. Despite the fact that the book is primarily intended for those 
who apply mathematics, the presentation is sufficiently careful to satisfy reasonable requirements as to 
rigor. 
The introductory chapter provides a brief but satisfactory background in the theory of functions of a 
real variable and the book continues through a detailed discussion of such topics as matrices, tensors, 
function theory and related subjects, asymptotic expansions, the classical equations (equations of heat 
conduction, wave equation, Laplace equation, etc.), and special functions. 

As in any book of this scope, there are a few notable omissions. The Fourier analysis might have been 
presented in more general form; the isoperimetric problems and direct methods of the variational calculus 
are absent, and, in the chapter on numerical methods, a wealth of detail on interpolation formuli leaves 
no room for many valuable topics, e.g. the Rayleigh-Ritz procedure and the techniques for finding eigen- 
values of matrices. It also seems to the reviewer that a more fortunate presentation of the Laplace Trans- 
form techniques might have been used. However, these latter remarks are greatly outweighed by the 
over-all value of the book which should find a place as an excellent reference for the mathematical physi- 


cist. 
G. F. CARRIER 
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Matrix and tensor calculus, with applications to mechanics, elasticity and aeronautics. 
By Aristotle D. Michal. Galcit Aeronautical Series. John Wiley & Sons, Inc. 
New York, 1947. xiii+132 pp. $3.00. 


The purpose of this book is to furnish engineering students (especially professional engineers) with 
an introduction to the theory and applications of the matrix and tensor analysis. In view of this purpose 
and the vast literature in each of these subjects, the author develops only: (1) those aspects of the theory 
of matrices which are of immediate application to the problems of vibration theory; (2) those elements 
of tensor analysis which suffice for a tensor formulation of the Lagrange, elasticity and hydrodynamical 
equations. The mathematical theory is presented in concise form and some applications of the theory 
to practice are clearly indicated. The problems supplement the material of the text and, though few in 
number, should provide the engineer with some of the necessary practice in mastering the new tools. 

Some indication of the material contained in the book may be obtained by a brief summary of the 
topics covered. Part I deals with such topics as addition, multiplication, differentiation and integration 
of matrices, the Cayley-Hamilton theorem as applied to the calculation of the inverse of a given matrix, 
systems of matrix differential equations and the solution of oscillation problems for characteristic fre- 
quencies and amplitudes. The problems contain some applications to flutter theory. Most of Part II is 
concerned with tensor analysis in Euclidean three space. Scalers, vectors, tensors (contravariant and co- 
variant) are defined. As examples of tensors, the author studies the metric and strain tensors. The La- 
place, Poisson, and wave equations are derived in general curvilinear coordinates. The related problem 
of the connection between the modern and classical theories of tensors in orthogonal curvilinear co- 
ordinates (see L. Brillouin, Les Tenseurs, Dover, New York, 1946, p. 114) is not discussed in the text. 
However, the author furnishes a detailed and interesting tensor treatment of elasticity theory. The 
last two chapters of the book deal with the tensor analysis of n-dimensional Riemann space. Geodesics 
and the Lagrange equations of motion are discussed. Finally, the Navier-Stokes equations for viscous 
fluids are formulated in tensor form and a treatment of the boundary layer problem by a perturbation 
method is indicated (the author gives credit to Dr. C. C. Lin for this latter work). 

In view of the many uses of matrix theory and the value of tensor analysis in gaining an insight 
into the geometric structure of partial differential equations, it has become essential that the engineer 
obtain some grasp of these subjects. This text should furnish the engineer with a good introduction to these 
fields. 

N. CoBuRN 


Introduction to the theory of equations. Second edition. By Lois W. Griffiths, John 
Wiley and Sons, Inc., New York; Chapman and Hall, Ltd., London, 1947. 278 
pp. $3.50 


The usual topics are treated here with unusual attention to detail. Each general theorem is pre- 
ceded by illustrative special cases, and numerous problems are included. Proofs of the main theorems 
are accomplished by gradually building up the desired statements out of simpler ones. The author is as 
fully concerned with teaching the student how to prove these theorems as with the theorems themselves. 
The result is a somewhat long but very precise text which might best serve the student who intends to 
continue in mathematics. 

H. J. GREENBERG 


Statistical analysis in biology. Second Edition. By K. Mather. Interscience Publishers, 
Inc., New York, 1947. 267 pp. $5.00. 


This book is written for biologists with a rather weak background in mathematics. This aim has been 
so well carried out that many mathematicians will find the lack of analytical short-cuts annoying. They 
have no ground for complaint. The book can be recommended to biologists as one of a small number 
with modern techniques and relatively little mathematics. It can be recommended to mathematicians 
interested in how statistics is actually used—in biology at large, and more especially, in the author's 
field of genetics. 

The second edition differs from the first by the correction of some errors and the addition of a thir- 
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teenth chapter on the angular and probit transformations. The book’s special features include (7) the only 
available careful, elementary account of single degrees of freedom and (ii) the use of small sample tests 
before their large sample brethren. The book seems remarkably free of typographical errors, the re- 
viewer noting only “ten” for “nine” near the top of pages 228 and 231. 

Although it is a revised edition, there are ten or a dozen points with which the reviewer differs. Some 
are errors of fact, others are differences in philosophy. Some of the most important seem to have arisen 
by incautious reading of R. A. Fisher—paying heed to what Fisher’s words seem to mean rather than to 
what Fisher was very careful to say. 

The “rigour” of statistical inferences is emphasized at page 11, without any warning that the exact- 
ness of these Fisherian procedures (which the reviewer uses himself) depends or the exact normality of 
the distributions found in practice—something difficult if not impossible to verify by experiment. 

A minor point on page 35 concerns kurtosis. Here overprevalence in the center is given equal 
weight with overprevalence in the tails as a cause of kurtosis. This is only half the classical error, which 
emphasized the center even more. 

The reader is advised on page 45 to put the larger mean square in the denominator of F. This is not 
consistent with the modern attitude toward the analysis of variance (e.g. Eisenhart, Biometrics 3 (1947), 
pp. 1-21). ee 

On page 58 the reader has completed the first analysis of the tomato data. He has found that A 
“must be supposed to... (be) greater than 1,” that B “does not... differ significantly from 1,” and 
“that there is no reason to suppose that a real difference . . . exists (between A and B).” (Parenthetic 
words ours.) The reviewer feels that the reader deserves to have the situation clarified. 

The reader who reads on page 76 “is... appropriate ...in this case, it being usual to take the 
highest order interaction for this purpose, unless . . . ” will be surprised to learn nine pages later, on page 
85, that he should use a different error term. 

The reviewer cannot combine the comments on pages 78 and 131 on pooling estimates of error into 
a coherent philosophy. His own philosophy is not coherent either, but he would do the opposite in both 
cases. Perhaps Mather’s third edition will clarify this problem for all of us. 

The trout fly example on page 123 seems somewhat incautiously chosen—the distributions seem 
so far from normality as to make the author’s concern about differences between significance at the 
4%, 2% or 1% level somewhat exaggerated. The data thus do not make an important point with the 


definiteness the point deserves. 

A major error of fact occurs on page 205, and again on page 212, where Mather states that maximum 
likelihood estimates have variances equal to 1/ni in samples of m, The example of the Cauchy distribution 
for n=1 or n=2 shows the falsity of this statement—here 7 is finite but the variance is infinite. It is known 
that as n—> © the variance is asymptotically 1/ni, but this is all. 

On page 213 we are told that “It will be readily appreciated that inconsistent statistics are utterly 


misleading and should not be used under any circumstances.” This seems exaggerated, in view of such 
inconsistent statistics for the mean of a normal distribution as the following: “Take the mean of the 
sample, or, if there are more observations than electrons in the metagalaxy, take twice the mean.” (This 


is essentially due to Mosteller.) 

On page 222 we find turned into n—1 ina most ad hoc way. It is regrettably true that the maximum 
likelihood estimate of the variance of a normal population involves division by n, and not n—1. Regret- 
table, since this encourages n-dividers in their recreant ways. But since it is so, it should be admitted. 

On page 239 we are told that the variance of #/?, x Poisson, is exactly independent of the mean. 
Asymptotic and exact statements are confounded here, as are the means of the population and sample. 

Joun W. TUKEY 

















